Skip to main content

Thank you for visiting You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Magic in twisted transition metal dichalcogenide bilayers


The long-wavelength moiré superlattices in twisted 2D structures have emerged as a highly tunable platform for strongly correlated electron physics. We study the moiré bands in twisted transition metal dichalcogenide homobilayers, focusing on WSe2, at small twist angles using a combination of first principles density functional theory, continuum modeling, and Hartree-Fock approximation. We reveal the rich physics at small twist angles θ < 4, and identify a particular magic angle at which the top valence moiré band achieves almost perfect flatness. In the vicinity of this magic angle, we predict the realization of a generalized Kane-Mele model with a topological flat band, interaction-driven Haldane insulator, and Mott insulators at the filling of one hole per moiré unit cell. The combination of flat dispersion and uniformity of Berry curvature near the magic angle holds promise for realizing fractional quantum anomalous Hall effect at fractional filling. We also identify twist angles favorable for quantum spin Hall insulators and interaction-induced quantum anomalous Hall insulators at other integer fillings.


In condensed matter physics, simple and elegant models have often brought new ideas and started new paradigms. Celebrated examples include the Hubbard model for strongly correlated electron system1, the Tomonaga-Luttinger model for one-dimensional electron liquid2,3, and the Kitaev model for non-Abelian quantum spin liquid4, to name a few. As toy models are designed to illustrate key concepts in the simplest form, they are rarely realized directly in real materials, whose atomic-scale electronic structures are inevitably more complex. The recent advent of long-wavelength moiré superlattices based on 2D van der Waals structures provides a new and promising venue for physical realization and quantum simulation of model Hamiltonians. In magic-angle twisted bilayer graphene5 (TBG), experiments have discovered a variety of correlated electron states6,7,8,9,10,11,12 facilitated by flat moiré bands.

More recently, moiré superlattices of semiconducting transition metal dichalcogenides (TMD) have attracted interest as a potentially simpler and more robust platform for simulating the Hubbard model on an emergent lattice13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32. Each lattice site represents a low-energy electronic orbital in the moiré unit cell that spreads over many atoms. These semiconductor moiré systems can thus be viewed as artificial 2D solids—a periodic array of “magnified atoms”19. The atomic potential depth and interatomic bonding are highly tunable by the choice of TMD materials, the twist angle13 and the displacement field20,23. Thus, TMD-based moiré materials provide a favorable platform for simulating idealized models in two dimensions.

In this work, we predict the realization of generalized Kane−Mele models with topological flat band, interaction-driven Haldane insulator and Mott insulators in twisted TMD homobilayers at small twist angles. Contrary to current thoughts, we show by band structure calculation and analytical derivation that a magic twist angle exists in twisted TMD homobilayers, where the topmost valence miniband from the ±K-valleys is almost perfectly flat and well separated from other bands. This band carries a spin/valley Chern number and is well described by a generalized Kane−Mele model33.

At half filling of this topological flat band, we show that repulsive interactions drives spontaneous spin/valley polarization leading to Haldane’s quantum anomalous Hall insulator34. We further find an out-of-plane displacement field drives a transition from the Haldane insulator into a Mott insulator. Depending on the twist angle, this Mott state is either a spin/valley polarized ferromagnet or features intervalley coherence that spontaneously breaks the spin/valley U(1) symmetry. Thus our work reveals a rich phase diagram of topological, correlated and broken-symmetry insulators enabled by the flat band in TMD homobilayers at small twist angles below the 4–5 range in current experimental studies23,35.

Due to spin-valley locking36, monolayer TMDs such as WSe2 and MoTe2 feature top valence bands with spin- at +K valley and spin- at −K. We study TMD homobilayers with a small twist angle θ starting from AA stacking, where every metal (M) or chalcogen (X) atom on the top layer is aligned with the same type of atom on the bottom layer. In such twisted structure, the K points of the two layers are slightly displaced and form the two corners of the moiré Brillouin zone, denoted as κ±. A set of spin- () moiré bands is formed from hybridized +K (−K) valley bands of the two layers. The complete filling of a single moiré band including spin degeneracy thus requires 2 holes per moiré unit cell.


In order to obtain accurate moiré band structures, we perform large-scale density functional theory calculations with the SCAN+rVV10 van der Waals density functional37, which captures the intermediate-range vdW interaction through its semilocal exchange term. Focusing on twisted bilayer WSe2, we find that lattice relaxation has a dramatic effect on moiré bands. Our DFT calculations at θ = 5.08 with 762 atoms per unit cell show a significant variation of the layer distance d in different regions on the moiré superlattice, as shown in Fig. 1b. d = 6.7 Å is smallest in MX and XM stacking regions, where the metal atom on top layer is aligned with chalcogen atom on the bottom layer and vice versa, while d = 7.1 Å is largest in MM region where metal atoms of both layers are aligned. With the fully relaxed structure, the low-energy moiré valence bands of twisted bilayer WSe2 are found to come from the ±K valley (shown in Fig. 1c), as opposed to the Γ valley in previous computational studies38 and consistent with recent works35,39,40.

Fig. 1: Comparison with large-scale DFT calculations.
figure 1

a The κ± points of the moiré Brillouin zone are formed from the K points of the monolayer Brillouin zones, which are rotated by ±θ/2. b The interlayer distance of the twisted WSe2 structure obtained from DFT is shown, demonstrating a large variation between the MM and XM/MX regions. c The continuum band structure (blue lines) is plotted in comparison with large-scale DFT calculations (black dots) at twist angle θ = 5.08, showing excellent agreement.

At small twist angles, the large size of moiré unit cell makes it difficult to perform DFT calculations directly on twisted TMD homobilayers. An alternative and complementary approach, introduced by Wu et al.14, is the continuum model based on an effective mass description, which models the formation of moiré bands using spatially modulated interlayer tunneling ΔT(r) and layer-dependent potential Δ1,2(r). The continuum model Hamiltonian for ±K valley bands is given by

$${{{{{{{{\mathcal{H}}}}}}}}}_{\uparrow }=\left(\begin{array}{ll}-\frac{{\hslash }^{2}{({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\boldsymbol{\kappa }}}}}}}}}_{+})}^{2}}{2{m}^{* }}+{{{\Delta }}}_{1}({{{{{{{\bf{r}}}}}}}})&{{{\Delta }}}_{\rm {T}}({{{{{{{\bf{r}}}}}}}})\\ {{{\Delta }}}_{\rm {T}}^{{{{\dagger}}} }({{{{{{{\bf{r}}}}}}}})&-\frac{{\hslash }^{2}{({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\boldsymbol{\kappa }}}}}}}}}_{-})}^{2}}{2{m}^{* }}+{{{\Delta }}}_{2}({{{{{{{\bf{r}}}}}}}})\end{array}\right)$$

and \({{{{{{{{\mathcal{H}}}}}}}}}_{\downarrow }\) as its time-reversal conjugate.

The continuum model approach is valid at small twist angle where the moiré wavelength is large enough. In this case, the atom configuration within any local region of a twisted bilayer is identical to that of an untwisted bilayer with one layer laterally shifted relative to the other by a corresponding displacement vector d0. For example, \({{{{{{{{\bf{d}}}}}}}}}_{0}=0,-\left({{{{{{{{\bf{a}}}}}}}}}_{1}+{{{{{{{{\bf{a}}}}}}}}}_{2}\right)/3,\left({{{{{{{{\bf{a}}}}}}}}}_{1}+{{{{{{{{\bf{a}}}}}}}}}_{2}\right)/3\), with a1,2 the primitive lattice vector of a monolayer, correspond to the MM, MX and XM regions respectively. Therefore, the moiré potentials for twisted TMD bilayers ΔT(r) and Δ1,2(r) as a function of coordinate on the moiré superlattice can be determined from the valence band edges of the untwisted bilayer as a function of the corresponding shift vector20. In the lowest harmonic approximation, ΔT(r) and Δ1,2(r) are sinusoids that interpolate between MM, MX and XM regions14:

$${{{\Delta }}}_{1,2}({{{{{{{\bf{r}}}}}}}})=2V\mathop{\sum}\limits_{j=1,3,5}\cos ({{{{{{{{\bf{g}}}}}}}}}_{j}\cdot {{{{{{{\bf{r}}}}}}}}\pm \psi )$$
$${{{\Delta }}}_{\rm {T}}({{{{{{{\bf{r}}}}}}}})=w(1+{e}^{-i{{{{{{{{\bf{g}}}}}}}}}_{2}\cdot {{{{{{{\bf{r}}}}}}}}}+{e}^{-i{{{{{{{{\bf{g}}}}}}}}}_{3}\cdot {{{{{{{\bf{r}}}}}}}}})$$

where gj are (j − 1)π/3 counter-clockwise rotations of the moiré reciprocal lattice vector \({{{{{{{{\bf{g}}}}}}}}}_{1}=(4\pi \theta /\sqrt{3}{a}_{0},0)\), and a0 is the monolayer lattice constant. Up to an overall energy scale, the continuum model depends only on the dimensionless parameters \(\alpha \equiv V{\theta }^{2}/({m}^{* }{a}_{0}^{2})\), w/V and ψ.

From our DFT calculation for untwisted bilayers with relaxed layer distance and using the vacuum level as an absolute reference energy for the band edge, we obtain the continuum model parameters V = 9.0 meV, ψ = 128, and w = 18 meV. Importantly, the interlayer tunneling strength w is twice larger than previously reported14. To demonstrate the accuracy of the continuum model method, we compare in Fig. 1c the band structures computed by large-scale DFT directly at θ = 5.08 and by the continuum model with the above parameters, finding excellent agreement. Details on the DFT calculation can be found in Supplementary Note 1 (see supplemental material for details on DFT calculations and continuum model).

We remark that different approaches31,32,40 can lead to different conclusions on topology. Thus, we utilize a method to determine band topology directly from our large-scale DFT band structure based on symmetry eigenvalues. As detailed in Supplementary Note 2 (see supplemental material for details of topology of DFT wavefunctions), we are able to isolate bands from the ±K valley and compute their C3z eigenvalue at the high symmetry momenta γ, κ±, which determine their Chern number (mod 3)41. The C3z eigenvalues for the first two bands, summarized in Table 1, are consistent with the first two bands having non-trivial valley Chern number \({{{{{{{{\mathcal{C}}}}}}}}}_{K,1}={{{{{{{{\mathcal{C}}}}}}}}}_{K,2}=1\).

Table 1 C3z eigenvalues of the first two bands from each valley, computed from large-scale DFT wavefunctions at high symmetry momentum points.

Using the new continuum model parameters established above, along with the lattice constant a0 = 3.317 Å42 and the effective mass m* = 0.43me43,44, we calculate the band structure of twisted bilayer WSe2, Ei(k), at various twist angles, as shown in Fig. 2a. The bandwidth of the first band, \(W=\mathop{\max }\nolimits_{{{{{{{{\bf{k}}}}}}}}}{E}_{1}({{{{{{{\bf{k}}}}}}}})-\mathop{\min }\nolimits_{{{{{{{{\bf{k}}}}}}}}}{E}_{1}({{{{{{{\bf{k}}}}}}}})\), as well as the (direct or indirect) band gaps εij between pairs of bands (i, j) = (1, 2) and (2, 3), \({\varepsilon }_{ij}=\mathop{\min }\nolimits_{{{{{{{{\bf{k}}}}}}}}}{E}_{i}({{{{{{{\bf{k}}}}}}}})-\mathop{\max }\nolimits_{{{{{{{{\bf{k}}}}}}}}}{E}_{j}({{{{{{{\bf{k}}}}}}}})\), is shown in Fig. 3. Focusing on topological features of the first two valence bands, we can divide the moiré band structure into three main regimes divided by θ1 ≈ 1.5 and θ2 ≈ 3.3.

Fig. 2: Continuum model band structure and Berry curvature at various twist angles.
figure 2

a The band structure Ei(k) along with the Chern numbers of the first two bands and (b) the (scaled) Berry curvature \(| {\kappa }_{+}{| }^{2}{{{{{{{\mathcal{F}}}}}}}}({{{{{{{\bf{k}}}}}}}})\) of the first band is shown for the continuum model at θ = 1, 1.43, 1.67, and 2.5. At θ = 2.5, the first band maxima is located at the κ± points and the Berry curvature is peaked at γ. At θ = 1, the band maximum is instead at γ and \({{{{{{{\mathcal{F}}}}}}}}\) is peaked at κ±. During the crossover region between these angles, E1(k) and \({{{{{{{\mathcal{F}}}}}}}}({{{{{{{\bf{k}}}}}}}})\) both become very flat. We find that the band dispersion E1(k) is flattest at θ ≈ 1.43 and the Berry curvature \({{{{{{{\mathcal{F}}}}}}}}({{{{{{{\bf{k}}}}}}}})\) is most uniform at θ ≈ 1.67, both shown.

Fig. 3: The bandwidth of the first band W, and indirect band gap between the first two pairs of bands ε12 and ε23.
figure 3

The bandwidth is minimized at θm, while being well separated from the remaining bands. The Chern numbers of the first two bands, \([{{{{{{{{\mathcal{C}}}}}}}}}_{K,1},{{{{{{{{\mathcal{C}}}}}}}}}_{K,2}]\), is shown before and after the ε23 gap closing at θ = θ1. For θ ≥ θ2, ε12 vanishes.

First, for θ < θ1, the top two bands are well separated from the rest of the spectrum, and carry opposite Chern number \([{{{{{{{{\mathcal{C}}}}}}}}}_{K,1},{{{{{{{{\mathcal{C}}}}}}}}}_{K,2}]=[+1,-1]\). The bandwidth of the first band W < 1 meV remains very small throughout. In this regime of very small twist angles, the character of the top two valence bands can be understood from an effective tight binding model on a moiré honeycomb lattice that takes the form of a Kane−Mele model, as suggested in the insightful work of Wu et al.14 As we will later show, the original Kane−Mele description with up to second nearest neighbor hopping terms only describes the band structure well for very small angles θ 1. As θ increases towards θ1, longer range hopping terms become more important.

At θ = θ1, the band gap ε23 closes and the Chern number of the top two bands changes to [+1, +1]. In this second regime, θ1 < θ < θ2, both top bands have same Chern number [+1, +1] and are still all separated by a sizable gap ε12, ε23 > 0. The bandwidth of the first band increases rapidly with θ, reaching around W ≈ 20 meV at θ2 (not shown). Finally, in the third regime, θ2 < θ 5.4, the indirect gap ε12 vanishes, but the direct gap remains open. The Chern number of the top two bands remains well defined at [+1, +1], but the bands now overlap in energy and are highly dispersive. In both the second and third regimes, ε23 > 0; thus, the top two bands together form a gapped \({{{{{{{\mathcal{C}}}}}}}}=2\) manifold. Beyond θ 5.4, the gap ε23 also vanishes and the top two bands are no longer isolated (Supplementary Note 3; see supplemental material for derivation of the analytic magic angle). Topology of the continuum model at θ ≈ 5 is consistent with that determined directly from large-scale DFT in Table 1, further strengthening our confidence in the continuum model description even up to larger angles.

For θ < θ2, especially near θ2 where the first band is more dispersive, the spin Chern number \({{{{{{{\mathcal{C}}}}}}}}=1\) and ε12 > 0 is favorable for a quantum spin Hall insulator at a filling of n = 2 holes per moiré unit cell. Also, for the wide range of angles θ1 < θ 5.4, ε23 > 0 and the top two bands both carry spin Chern number \({{{{{{{\mathcal{C}}}}}}}}=1\), giving rise to a double quantum spin Hall state with two sets of counter-propagating spin-polarized edge modes at n = 4.

We now address the bandwidth W, which shows a sharp minimum at θ = θm ≈ 1.43 reminiscent of the magic angle in TBG. To understand this, notice that the top band, shown in Fig. 2a, has two qualitatively different behaviors in the small and large θ limit. For θ 2.5, E1(k) reaches its maximum at κ± and minimum at γ, which can be understood from the weak moiré effects at small α. For small θ 1, the opposite holds and E1(k) is maximum is at γ and minimum at κ±, which can be understood from the effective Kane−Mele model, which we will derive explicitly. At the crossover between these two limits, the band maxima and minima must switch locations in the moiré Brillouin zone, potentially leading to a flat band. As can be clearly seen, W achieves a minimum at a particular magic angle θm during this crossover. At θm the gap to the next state ε12 ≈ 3.7 meV is much larger than the bandwidth W ≈ 0.1 meV. The band structure at θm is shown in Fig. 2a, which shows that the first band is almost completely flat and separated from the next band. For even smaller θ, both ε12 and W vanish, but the ratio W/ε12 diverges. Thus, we may view θm as the angle at which the top band is flattest while still being well isolated from the rest of the spectrum.

Analytic progress can be made in estimating the magic angle by considering the dispersion near γ. Assuming that the bandwidth will be minimized near the angle at which E1(γ) changes from minima to maxima, expanding \({E}_{1}({{{{{{{\boldsymbol{\gamma }}}}}}}}+{{{{{{{\bf{k}}}}}}}})\approx {E}_{1}({{{{{{{\boldsymbol{\gamma }}}}}}}})+\frac{{{{{{{{{\bf{k}}}}}}}}}^{2}}{2{m}_{\gamma }}+{{{{{{{\mathcal{O}}}}}}}}({{{{{{{{\bf{k}}}}}}}}}^{3})\), the effective mass mγ should diverge near the crossover. Let \({\tilde{\theta }}_{m}\) to be the angle at which \({m}_{\gamma }^{-1}=0\). Then, considering only the six most relevant states at γ, we have (Supplementary Note 4; see supplemental material for continuum model band structure at higher twist angles)

$${\tilde{\theta }}_{m}^{-2}=\frac{8{\pi }^{2}}{9{m}^{* }{a}_{0}^{2}}\left(\frac{1}{{{{{{{{{\mathcal{E}}}}}}}}}_{{n}_{0}}-{{{{{{{{\mathcal{E}}}}}}}}}_{{n}_{0}+1}}+\frac{1}{{{{{{{{{\mathcal{E}}}}}}}}}_{{n}_{0}}-{{{{{{{{\mathcal{E}}}}}}}}}_{{n}_{0}-1}}\right)$$

where \({{{{{{{{\mathcal{E}}}}}}}}}_{n}=2w\cos (\pi n/3)+2V\cos (2\pi n/3-\psi )\), and n0 is the integer (mod 6) which maximizes \({{{{{{{{\mathcal{E}}}}}}}}}_{n}\) (n0 = 1 for WSe2 parameters). We find that Eq. (4) provides a decent estimate for θm in the cases considered. In WSe2, we have \({\tilde{\theta }}_{m}=1.4{7}^{\circ }\), compared to θm = 1.43 at which the bandwidth is minimized.

Next, we turn to the Berry curvature \({{{{{{{\mathcal{F}}}}}}}}({{{{{{{\bf{k}}}}}}}})\) of the top band, shown in Fig. 2b. In all cases, the first band has Chern number \({{{{{{{{\mathcal{C}}}}}}}}}_{K,1}=\frac{1}{2\pi }{\int}_{{{{{{{{\rm{BZ}}}}}}}}}{{{{{{{\mathcal{F}}}}}}}}{\rm {d}}{{{{{{{\bf{k}}}}}}}}=1\); however, the distribution changes very drastically as θ is varied. At θ 2, \({{{{{{{\mathcal{F}}}}}}}}\) is peaked around the band minimum at γ. At θ 1, \({{{{{{{\mathcal{F}}}}}}}}\) is instead sharply peaked at the κ± points. Near the crossover region, the distribution of \({{{{{{{\mathcal{F}}}}}}}}\) shifts from γ to κ±, and can become very evenly distributed. We find that \({{{{{{{\mathcal{F}}}}}}}}\) is most evenly distributed near θ = 1.67, shown in Fig. 2b, where \({{{{{{{\mathcal{F}}}}}}}}\) becomes almost uniform in the Brillouin zone. The uniform distribution of \({{{{{{{\mathcal{F}}}}}}}}\) is reminiscent to that of Landau levels. Time-reversal symmetry forces the corresponding spin- bands from the −K valley to have opposite Chern number.

We emphasize that the physics of the magic angle arises due to the crossover between two qualitatively different behaviors of the first band at low and high angles. Additional factors unaccounted for by the continuum model may result in, for example, angle-dependent model parameters. However, as long as the qualitative behaviors at small and large angles are unchanged, there will still be crossover regime at which the band becomes flat. Even when the bands are not perfectly flat, a diverging mass can still give rise to a diverging higher-order van Hove singularity23.

Recall that for θ < θ1, the top two bands carry opposite Chern number and are separated from the rest of the spectrum, suggesting a description in terms of an effective tight binding model. We now focus on θ < θ1 and directly derive an effective tight binding model for the first two moiré bands by explicitly constructing a basis of localized Wannier states. These Wannier states are constructed via a simple procedure which manifestly preserves the symmetries of the twisted homobilayer. Given the single-particle eigenstates \(\{|{\phi }_{n,{{{{{{{\bf{k}}}}}}}}}\rangle \}\) of the continuum Hamiltonian (1) for each k in the mBZ, we first construct a superposition of the first two (n = 1, 2) eigenstates, \(|{\tilde{\phi }}_{n,{{{{{{{\bf{k}}}}}}}}}\rangle ={\sum }_{m = 1,2}{U}_{nm}^{({{{{{{{\bf{k}}}}}}}})}\left|{\phi }_{m{{{{{{{\bf{k}}}}}}}}}\right\rangle\) using a 2 × 2 unitary matrix \({U}_{nm}^{({{{{{{{\bf{k}}}}}}}})}\), which maximizes the layer polarization at every k:

$${P}_{{{{{{{{\bf{k}}}}}}}}}=\mathop{\sum}\limits_{n=1,2}{(-1)}^{n}\left\langle {\tilde{\phi }}_{n,{{{{{{{\bf{k}}}}}}}}}| ({{{{{{{{\mathcal{P}}}}}}}}}_{-}-{{{{{{{{\mathcal{P}}}}}}}}}_{+})| {\tilde{\phi }}_{n,{{{{{{{\bf{k}}}}}}}}}\right\rangle ,$$

where \({{{{{{{{\mathcal{P}}}}}}}}}_{\pm }\) is the projector on to the top/bottom layer, so that \(|{\tilde{\phi }}_{1,{{{{{{{\bf{k}}}}}}}}}\rangle\) is chosen to mostly consist of states in the top layer, and similarly for \(|{\tilde{\phi }}_{2,{{{{{{{\bf{k}}}}}}}}}\rangle\) on the bottom layer. This uniquely specifies \(|{\tilde{\phi }}_{n,{{{{{{{\bf{k}}}}}}}}}\rangle\) up to a phase, which we choose to be real and positive at the XM (n = 1) or MX (n = 2) stacking regions (Supplementary Note 5; see supplemental material for details on the derivation of localized Wannier functions and tight binding model). The Wannier states at moiré lattice vector R is then defined \(\left|{W}_{{{{{{{{\bf{R}}}}}}}}}^{n}\right\rangle =\frac{1}{\sqrt{{N}_{k}}}{\sum }_{{{{{{{{\bf{k}}}}}}}}}{e}^{-i{{{{{{{\bf{k}}}}}}}}\cdot {{{{{{{\bf{R}}}}}}}}}\left|{\tilde{\phi }}_{n{{{{{{{\bf{k}}}}}}}}}\right\rangle\). They are localized about their centers with a root-mean-square distance aW ≈ 5 nm, and are also mostly composed of states in one layer: \(\left\langle {W}_{{{{{{{{\bf{R}}}}}}}}}^{1}| {{{{{{{{\mathcal{P}}}}}}}}}_{+}| {W}_{{{{{{{{\bf{R}}}}}}}}}^{1}\right\rangle \approx 0.83\) is mostly in the top layer, and vice versa for \(\left|{W}_{{{{{{{{\bf{R}}}}}}}}}^{2}\right\rangle\).

It is straightforward to obtain the hopping matrix elements of the effective tight binding model in the Wannier basis for the top two bands as a function of θ. Figure 4b shows the nth nearest neighbor hopping matrix elements tn obtained in this way, up to n = 5. As anticipated, the effective tight binding model at θ < θ1, including the spin/valley degrees of freedom, is a generalized Kane−Mele model with sites centered on the honeycomb lattice formed by MX and XM stacking regions14.

Fig. 4: Wannier functions and tight binding model parameters.
figure 4

a Wannier functions at the magic angle, b tight binding parameters as a function of θ, and c the bandwidth of the top band in the effective tight binding model TBn keeping up to nth nearest neighbor hopping terms, compared to that of the continuum model.

The tight binding Hamiltonian is found to be

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{TB}}}}}}}}}={t}_{1}\mathop{\sum}\limits_{\langle i,j\rangle ,\sigma }{c}_{i\sigma }^{{{{\dagger}}} }{c}_{j\sigma }+| {t}_{2}| \mathop{\sum}\limits_{\langle \langle i,j\rangle \rangle ,\sigma }{e}^{i\phi \sigma {\nu }_{ij}}{c}_{i\sigma }^{{{{\dagger}}} }{c}_{j\sigma }+\cdots$$

where \({c}_{i\sigma }^{{{{\dagger}}} },{c}_{i\sigma }\) are fermionic creation/annihilation operators, σ = ± is the spin/valley degree of freedom, the sum 〈i, j〉 (〈〈i, j〉〉) is over (next) nearest neighboring sites i, j of the honeycomb lattice, and νij = ±1 depending on whether the path i → j turns right (+) or left (−). The parameter t1 is real, while t2 ≡ t2eiϕ is complex, and contain longer range hopping terms. We find that tn quickly reduce in magnitude with hopping distance n, and only the second neighbor hopping has a significant imaginary component. In Fig. 4c, we show the bandwidth of the top band, W, in the effective tight binding model TBn including up to tn hopping terms, compared to that of the continuum model. For θ 1, TB2 already captures the band structure very well. Near the magic angle, higher range hoppings become more important in capturing the flatness of the band.

For the small twist angles θ < θ1 considered, since the size of the Wannier orbitals are small compared to the moiré unit cell, the dominant interaction is a simple on-site Hubbard term \({{{{{{{{\mathcal{H}}}}}}}}}_{U}=U{\sum }_{i}{n}_{i\uparrow }{n}_{i\downarrow }\). We also estimate U ~ e2/(ϵaW) ≈ 70 meV at θm, using a realistic relative dielectric constant ϵ = 4 and aW =  5nm, which is significantly larger than the tight binding parameters tn. Therefore at such small twist angle, twisted WSe2 homobilayers are in the strong-coupling regime, in contrast with θ ~ 4–5 where the bandwidth is comparable to the interaction strength23,35.

In the following, we shall focus on the strongly correlated regime θ < θ1 at a filling of n = 1 holes per moiré unit cell, where we expect the flat bands will favor the quantum anomalous Hall (QAH) insulator due to spontaneous spin/valley polarization. The reason is as follows. First, spin/valley polarized states filling the top band of \({{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{TB}}}}}}}}}\) with σ = ± are exact eigenstates of our interacting model, because the spin-orbit coupling in \({{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{TB}}}}}}}}}\) conserves the z-spin/valley component. Then, these spin/valley polarized states avoid the Hubbard interaction, and also minimize the kinetic energy in the case of a completely flat top band. Minimizing both parts of the Hamiltonian, they necessarily are many body ground states of the model at n = 1 filling. The complete polarization of \({{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{TB}}}}}}}}}\) exactly corresponds to Haldane’s model for QAH insulator.

To include corrections coming from the finite bandwidth W of the flat band, we solve our interacting problem within the Hartree−Fock (HF) approximation, where the Hubbard interaction is decoupled as \({n}_{i\uparrow }{n}_{i\downarrow }\simeq \langle {n}_{i\uparrow }\rangle {n}_{i\downarrow }+\langle {n}_{i\downarrow }\rangle {n}_{i\uparrow }-\langle {c}_{i\uparrow }^{{{{\dagger}}} }{c}_{i\downarrow }\rangle {c}_{i\downarrow }^{{{{\dagger}}} }{c}_{i\uparrow }-\langle {c}_{i\downarrow }^{{{{\dagger}}} }{c}_{i\uparrow }\rangle {c}_{i\uparrow }^{{{{\dagger}}} }{c}_{i\downarrow }\), up to a constant term, and where the expectation values for the spin and density at each site are determined self-consistently by iteration. Our numerical solutions of the HF equations are shown in Fig. 5a as a function of twist angle and interaction strength. As expected, we observe a transition from a metallic state to a ferromagnetic QAH insulator polarized along z when U increases. Within HF, this transition can be understood as follows. The fully polarized states yield a rigid shift of the bands by σU/2. When U is larger than the non-interacting bandwidth W, a full gap opens and the ferromagnetic state fully fills one of Chern bands of \({{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{TB}}}}}}}}}\), which leads to a QAH phase. We remark that the appearance of the QAH phase relies on both the non-trivial Chern number as well as the fact that the band is flat and isolated, features which are maximized at the magic angle, as illustrated by a dip of the insulating phase above θ 1.4.

Fig. 5: Numerical solution of the self-consistent Hartree−Fock approximation.
figure 5

The phase diagram is shown as a function of the twist angle and (a) interaction strength at fixed displacement field Δ = 0, or b displacement field at fixed interaction strength U = 10t1, including up to fifth nearest neighbor hopping terms. The insulating phases denoted by QAH, FMz, 120 AFMxy, and FMxy are described in the main text. The dashed blue line shows the boundary of the FMz phases determined by the first magnon instability. The colors indicate the charge gap. Hartree−Fock calculations were done with a \(\sqrt{3}\times \sqrt{3}\) unit cell of 6 atom sites, and a 30 × 30 grid of (k) points.

To precisely locate the transition between QAH insulator and the metallic phase, we compute the magnon excitation spectrum above the fully ferromagnetic state by exact diagonalization (ED) of the interacting Hamiltonian projected on the spin-1 excitation subspace45. For large U, this spectrum is gapped and the QAH is robust against spin flips. Decreasing U eventually brings one magnon at zero energy, which destabilizes the ferromagnetic states and drives the transition to a metal. As shown in Fig. 5a, the ED results almost perfectly agrees with the HF boundaries, putting them on firmer grounds.

For the large values of U relevant to WSe2, the magnons have a large gap, and the lowest excitation corresponds to an interband transition between two bands with same spin. The QAH phase is thus protected by a gap ε12 ≈ 3.7 meV near the magic angle, leading to quantum Hall effect at elevated temperature.

We also highlight that the QAH may also be observed for larger twist angles, where the first band still carries a non-zero Chern number (Fig. 2), and its bandwidth remain small compared to the estimated U (Fig. 3). Likewise, the second band is topological and quite flat for θ ~ 2–3 (Supplementary Note 3), and therefore QAH may also be observed at a filling of n = 3. Twisted TMD bilayers with ±K-valley bands are thus expected to be an intrinsically robust platform for interaction-induced QAH phases.

It is interesting and worthwhile to compare the QAH phase in twisted TMD and graphene bilayers. Anomalous Hall effect and its quantization have been experimentally observed in magic-angle TBG11,12, where the alignment with hBN substrate is likely the origin of valley Chern number46 and both spin and valley degeneracy are lifted due to repulsive interaction in the flat band47,48,49,50,51. Due to the presence of SU(2)-invariant spin degrees of freedom, QAH in TBG is subject to the adverse effect of gapless thermal fluctuation, which forbids long-range order at finite temperature in the thermodynamic limit. In contrast, spin-valley locking in TMD systems enables robust Ising-type spin/valley order that leads to QAH effect at lower temperature.

Another great advantage of twisted TMD bilayer is their high degree of tunability, in particular with respect to applied electric fields19,35,52. Due to the layer polarization of the Wannier basis states, the displacement field can be modeled as a sublattice symmetry breaking term \({{{{{{{{\mathcal{H}}}}}}}}}_{{{\Delta }}}=\frac{{{\Delta }}}{2}{\sum }_{i}{s}_{i}{c}_{i\sigma }^{{{{\dagger}}} }{c}_{i\sigma }\), where si is (−)1 for i in the A (B) sublattice. Including this term in our HF treatment, we can investigate which phases should neighbor the QAH ferromagnet in experiments. We present our solutions of the HF equations as a function of twist angle and displacement field in Fig. 5b. There, we fix U = 10t1, a tradeoff between the large U of the homobilayer system and the convergence rate of the HF self-consistent iteration algorithm. We find it necessary to consider an enlarged \(\sqrt{3}\times \sqrt{3}\) unit cell in order to describe all ordered phases of the model.

At small displacement fields, the topmost moiré band remain relatively flat and our earlier arguments for spin/valley polarization still apply. This is confirmed by our HF solutions for Δ  2t1 (Fig. 5b), which exhibit full spin polarization along the z axis. In this region, a transition nevertheless occurs at \({{\Delta }}=6\sqrt{3}| {t}_{2}| \sin \phi\) (up to tn≥3 terms), where the single-particle gap between the two moiré bands closes, and their Chern numbers change from [+1, −1] to [0, 0]. This gap closing line marks the transition between a QAH insulator and a topologically trivial ferromagnet with spin/valley polarization (FMz)34. As the displacement field further increases, the bandwidth W also grows, which decreases the magnon gap (see discussion above). The spin/valley polarized phases eventually become unstable when the magnon gap closes, which can be seen with the very good agreement between the phase boundaries determined with HF and ED (Fig. 5b).

Beyond this spin-wave instability line, our results show the emergence of two new Mott insulating phases, where holes are mostly localized on the A sublattice, and their spin either form an antiferromagnetic pattern (120 AFMxy), or ferromagnetically align in the xy plane (FMxy). Their appearance is most easily understood for large displacement fields, where the physics becomes analogous to that of localized moments on the triangular A sublattice. In the regime t2t1 Δ, U relevant for our system, their coupling is described by an effective XXZ model with Dzyaloshinskii−Moriya (DM) interactions

$${{{{{{{{\mathcal{H}}}}}}}}}_{S}=\,\mathop{\sum}\limits_{{\langle i,j\rangle }_{B}}\,{J}_{\parallel }{s}_{i}^{z}{s}_{j}^{z}+{J}_{\perp }({s}_{i}^{x}{s}_{j}^{x}+{s}_{i}^{y}{s}_{j}^{y})+D\left[({{{{{{{{\bf{s}}}}}}}}}_{i}\times {{{{{{{{\bf{s}}}}}}}}}_{j})\cdot {{{{{{{\bf{z}}}}}}}}\right],$$

which is derived in Supplementary Note 6 (see supplemental material for the derivation of the effective spin Hamiltonian). The parameters of this effective spin model are given by

$${J}_{\parallel }=\frac{4| \tilde{t}{| }^{2}}{U}+{{{{{{{\rm{Re}}}}}}}}\left(\frac{4{t}_{1}^{2}{t}_{2}}{{{{\Delta }}}^{2}}\right),$$
$${J}_{\perp }={{{{{{{\rm{Re}}}}}}}}\left(\frac{4{\tilde{t}}^{2}}{U}+\frac{4{t}_{1}^{2}{t}_{2}}{{{{\Delta }}}^{2}}\right),$$
$$D={{{{{{{\rm{Im}}}}}}}}\left(\frac{4{\tilde{t}}^{2}}{U}+\frac{4{t}_{1}^{2}{t}_{2}}{{{{\Delta }}}^{2}}\right),$$

with \(\tilde{t}={t}_{2}+{t}_{1}^{2}/{{\Delta }}\). In Eq. 8, we have separated exchange terms coming from different physical processes. The first ones \(\propto {\tilde{t}}^{2}/U\) arise from nearest neighbor tunneling on the triangular A sublattice, while the others \(\propto {t}_{2}{({t}_{1}/{{\Delta }})}^{2}\) originate from loop-exchange on the honeycomb lattice that do not involve any double occupancy.

For twist angles θ 1, t2t1 and Eq. (7) reduces to an antiferromagnetic (AFM) Heisenberg model, where J = J > 0 are dominated by the nearest neighbor tunneling on the triangular lattice. This simplified triangular lattice description, valid for very small twist angles, has been proposed in earlier studies of Mott insulators in twisted TMDs24,26,53. It was shown to yield an antiferromagnetic phase that the small residual DM interaction pins in the xy plane. This is the origin of the AFMxy phase observed in Fig. 5b. We also note that the weak-coupling version of AFMxy phase—an intervalley-coherent \(\sqrt{3}\times \sqrt{3}\) density wave23—has been proposed for the correlated insulating state at n = 1 in twisted bilayer WSe2 at θ ~ 4–535.

For larger twist angles, t2 becomes substantial and we observe that J becomes negative for the realistic parameter U Δ, dominated by a third-order exchange process \(\propto {t}_{1}^{2}{t}_{2}\) on the honeycomb lattice without double occupancy. Then, the FMxy phase is favored as shown in Fig. 5b. The competition between AFMxy and FMxy phases can be analyzed by solving Eq. (7) for classical spins. This approach, detailed in the Supplementary Note 6 (see supplemental material for the derivation of the effective spin Hamiltonian), gives a transition between the two phases when \(| D| =-\sqrt{3}{J}_{\perp }\). For Δ = 5t1, this criterion yields a critical twist angle θ = 0.95, which roughly agrees with our HF results. We note that the ferromagnetic phase due to J < 0 does not appear in twisted TMDs based on simplified triangular lattice descriptions.

Finally, we comment on the effect of nearest neighbor repulsion Vi, jninj to our HF phase diagram. This term favors the layer polarized phases, such as the FMxy and AFMxy which appear at large Δ. Small V therefore narrows the range in Δ at which the QAH phase appears. For large V, there is a sharp transition at Δ = 0 between layer polarized Mott insulating phases, which can lead to the strong hysteretic behavior of Mott ferroelectricity20,54. The long-range component of interactions can be controlled by screening from nearby metallic layers. Multiple recent experiments55,56 on twisted WSe2 homobilayers in the presence of a nearby WSe2 monolayer report strong screening effects when the monolayer is doped. This raises the interesting possibility of a screening-induced transition between the QAH and Mott insulating phases.


Our phase diagram demonstrates the high experimental tunability of TMD twisted homobilayers, where the applied displacement field can tune between quantum anomalous Hall phase and Mott insulators involving three types of magnetic orders: FMz, FMxy, and AFMxy. Despite being electrically insulating, the xy-ordered Mott insulators support coherent magnon transport23, which can be detected by optical spin injection and spatial-temporal mapping recently developed for TMD bilayers57. The experimental feasibility of tuning and distinguishing between topologically different insulators at the same filling adds to the attractiveness and desirability of TMD-based moiré systems.

In parallel to our work on twisted TMD homobilayers, a breakthrough experiment led by Kin Fai Mak and Jie Shan discovered unexpectedly a QAH phase with spontaneous spin/valley polarization in a TMD heterobilayer MoTe2/WSe2 at the filling of n = 1 tuned by displacement field58,59. Large-scale DFT calculation and wavefunction analysis reveal two dispersive moiré bands forming the Kane−Mele model, suggestive of a similar origin of QAH as described here.

Looking forward, the remarkable flat Chern band we found, combined with the uniformity of Berry curvature, suggests that twisted TMD homobilayers near magic angle may be an ideal setting for observing a fractional quantum anomalous Hall state at zero magnetic field.

Data availability

The data needed to evaluate the conclusions in the paper are present in the paper and the Supplementary Material. The full dataset generated during this study, including relaxed lattice structure and band structure obtained from DFT, tight binding model parameters, and self-consistent HF solutions, has been deposited in the Zenodo database60. Additional data related to this paper are available from the corresponding authors upon reasonable request.


  1. 1.

    Hubbard, J. Electron correlations in narrow energy bands. Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 276, 238–257 (1933).

    ADS  Google Scholar 

  2. 2.

    Tomonaga, S.-I. Remarks on Bloch’s method of sound waves applied to many-Fermion problems. Prog. Theor. Phys. 5, 544–569 (1950).

    ADS  MathSciNet  Google Scholar 

  3. 3.

    Luttinger, J. M. An exactly soluble model of a many fermion system. J. Math. Phys. 4, 1154–1162 (1963).

    ADS  MathSciNet  CAS  Google Scholar 

  4. 4.

    Kitaev, A. Anyons in an exactly solved model and beyond. Ann. Phys. 321, 2–111 (2006).

    ADS  MathSciNet  MATH  CAS  Google Scholar 

  5. 5.

    Bistritzer, R. & MacDonald, A. H. Moiré bands in twisted double-layer graphene. Proc. Natl Acad. Sci. USA 108, 12233–12237 (2011).

    ADS  PubMed  PubMed Central  CAS  Google Scholar 

  6. 6.

    Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature 556, 80–84 (2018).

    ADS  PubMed  CAS  Google Scholar 

  7. 7.

    Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices. Nature 556, 43–50 (2018).

    ADS  PubMed  CAS  Google Scholar 

  8. 8.

    Yankowitz, M. et al. Tuning superconductivity in twisted bilayer graphene. Science 363, 1059–1064 (2019).

    ADS  PubMed  CAS  Google Scholar 

  9. 9.

    Lu, X. et al. Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene. Nature 574, 653–657 (2019).

    ADS  PubMed  CAS  Google Scholar 

  10. 10.

    Cao, Y. et al. Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene. Nature 583, 215–220 (2020).

    ADS  PubMed  CAS  Google Scholar 

  11. 11.

    Sharpe, A. L. et al. Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene. Science 365, 605–608 (2019).

    ADS  PubMed  CAS  Google Scholar 

  12. 12.

    Serlin, M. et al. Intrinsic quantized anomalous hall effect in a moiré heterostructure. Science 367, 900–903 (2020).

    ADS  PubMed  CAS  Google Scholar 

  13. 13.

    Wu, F., Lovorn, T., Tutuc, E. & MacDonald, A. H. Hubbard model physics in transition metal dichalcogenide moiré bands. Phys. Rev. Lett. 121, 026402 (2018).

    ADS  PubMed  CAS  Google Scholar 

  14. 14.

    Wu, F., Lovorn, T., Tutuc, E., Martin, I. & MacDonald, A. H. Topological insulators in twisted transition metal dichalcogenide homobilayers. Phys. Rev. Lett. 122, 086402 (2019).

    ADS  PubMed  CAS  Google Scholar 

  15. 15.

    Tang, Y. et al. Simulation of hubbard model physics in wse2/ws2 moiré superlattices. Nature 579, 353–358 (2020).

    ADS  Article  CAS  Google Scholar 

  16. 16.

    Regan, E. C. et al. Mott and generalized wigner crystal states in wse2/ws2 moiré superlattices. Nature 579, 359–363 (2020).

    ADS  Article  CAS  Google Scholar 

  17. 17.

    Shabani, S. et al. Deep moiré potentials in twisted transition metal dichalcogenide bilayers. Nat. Phys. 17, 720–725 (2021).

    CAS  Google Scholar 

  18. 18.

    Jin, C. et al. Stripe phases in wse2/ws2 moiré superlattices. Nat. Mater. 20, 940–944 (2021).

    ADS  PubMed  CAS  Google Scholar 

  19. 19.

    Zhang, Y., Yuan, N. F. Q. & Fu, L. Moiré quantum chemistry: charge transfer in transition metal dichalcogenide superlattices. Phys. Rev. B 102, 201115 (2020).

    ADS  CAS  Google Scholar 

  20. 20.

    Zhang, Y., Liu, T. & Fu, L. Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers. Phys. Rev. B 103, 155142 (2021).

    ADS  CAS  Google Scholar 

  21. 21.

    Slagle, K. & Fu, L. Charge transfer excitations, pair density waves, and superconductivity in moiré materials. Phys. Rev. B 102, 235423 (2020).

    ADS  CAS  Google Scholar 

  22. 22.

    Xu, Y. et al. Correlated insulating states at fractional fillings of moiré superlattices. Nature 587, 214–218 (2020).

    ADS  PubMed  CAS  Google Scholar 

  23. 23.

    Bi, Z. & Fu, L. Excitonic density wave and spin-valley superfluid in bilayer transition metal dichalcogenide. Nat. Commun. 12, 642 (2021).

    PubMed  PubMed Central  CAS  Google Scholar 

  24. 24.

    Pan, H., Wu, F. & Das Sarma, S. Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer wse2. Phys. Rev. Res. 2, 033087 (2020).

    CAS  Google Scholar 

  25. 25.

    Morales-Durán, N., Potasz, P. & MacDonald, A. H. Metal-insulator transition in transition metal dichalcogenide heterobilayer moiré superlattices. Phys. Rev. B 103, 241110 (2021).

  26. 26.

    Zang, J., Wang, J., Cano, J. & Millis, A. J. Hartree-fock study of the moiré hubbard model for twisted bilayer transition metal dichalcogenides. Phys. Rev. B 104, 075150 (2021).

  27. 27.

    Padhi, B., Chitra, R. & Phillips, P. W. Generalized wigner crystallization in moiré materials. Phys. Rev. B 103, 125146 (2021).

    ADS  CAS  Google Scholar 

  28. 28.

    Zhai, D. & Yao, W. Theory of tunable flux lattices in the homobilayer moiré of twisted and uniformly strained transition metal dichalcogenides. Phys. Rev. Mater. 4, 094002 (2020).

    CAS  Google Scholar 

  29. 29.

    Cazalilla, M. A., Ochoa, H. & Guinea, F. Quantum spin hall effect in two-dimensional crystals of transition-metal dichalcogenides. Phys. Rev. Lett. 113, 077201 (2014).

    ADS  PubMed  CAS  Google Scholar 

  30. 30.

    Zhang, Z. et al. Flat bands in twisted bilayer transition metal dichalcogenides. Nat. Phys. 16, 1093–1096 (2020).

    CAS  Google Scholar 

  31. 31.

    Magorrian, S. J. et al. Multifaceted moiré superlattice physics in twisted wse2 bilayers. Phys. Rev. B 104, 125440 (2021).

  32. 32.

    Tang, H., Carr, S. & Kaxiras, E. Geometric origins of topological insulation in twisted layered semiconductors. Phys. Rev. B 104, 155415 (2021).

  33. 33.

    Kane, C. L. & Mele, E. J. Quantum spin hall effect in graphene. Phys. Rev. Lett. 95, 226801 (2005).

    ADS  PubMed  CAS  Google Scholar 

  34. 34.

    Haldane, F. D. M. Model for a quantum hall effect without Landau levels: condensed-matter realization of the “parity anomaly”. Phys. Rev. Lett. 61, 2015–2018 (1988).

    ADS  PubMed  CAS  Google Scholar 

  35. 35.

    Wang, L. et al. Correlated electronic phases in twisted bilayer transition metal dichalcogenides. Nat. Mater. 19, 861–866 (2020).

    PubMed  CAS  Google Scholar 

  36. 36.

    Xiao, D., Liu, G.-B., Feng, W., Xu, X. & Yao, W. Coupled spin and valley physics in monolayers of mos2 and other group-vi dichalcogenides. Phys. Rev. Lett. 108, 196802 (2012).

    ADS  PubMed  Google Scholar 

  37. 37.

    Peng, H., Yang, Z.-H., Perdew, J. P. & Sun, J. Versatile van der waals density functional based on a meta-generalized gradient approximation. Phys. Rev. X 6, 041005 (2016).

    Google Scholar 

  38. 38.

    Naik, M. H. & Jain, M. Ultraflatbands and shear solitons in moiré patterns of twisted bilayer transition metal dichalcogenides. Phys. Rev. Lett. 121, 266401 (2018).

    ADS  PubMed  CAS  Google Scholar 

  39. 39.

    Vitale, V., Atalar, K., Mostofi, A. A. & Lischner, J. Flat band properties of twisted transition metal dichalcogenide homo-and heterobilayers of mos _2, mose _2, ws _2 and wse _2. 2D Materials 8, 045010 (2021).

  40. 40.

    Kundu, S., Naik, M. H., Krishnamurthy, H.R. & Jain, M. Flat bands in twisted bilayer wse _2 with strong spin−orbit interaction. Preprint at (2021).

  41. 41.

    Fang, C., Gilbert, M. J. & Bernevig, B. A. Bulk topological invariants in noninteracting point group symmetric insulators. Phys. Rev. B 86, 115112 (2012).

    ADS  Google Scholar 

  42. 42.

    Mounet, N. et al. Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds. Nat. Nanotechnol. 13, 246–252 (2018).

    ADS  PubMed  CAS  Google Scholar 

  43. 43.

    Fallahazad, B. et al. Shubnikov–de haas oscillations of high-mobility holes in monolayer and bilayer wse 2: Landau level degeneracy, effective mass, and negative compressibility. Phys. Rev. Lett. 116, 086601 (2016).

    ADS  PubMed  Google Scholar 

  44. 44.

    Rasmussen, F. A. & Thygesen, K. S. Computational 2d materials database: electronic structure of transition-metal dichalcogenides and oxides. J. Phys. Chem. C 119, 13169–13183 (2015).

    CAS  Google Scholar 

  45. 45.

    Su, X.-F., Gu, Z.-L., Dong, Z.-Y. & Li, J.-X. Topological magnons in a one-dimensional itinerant flatband ferromagnet. Phys. Rev. B 97, 245111 (2018).

    ADS  CAS  Google Scholar 

  46. 46.

    Song, J. C. W., Samutpraphoot, P. & Levitov, L. S. Topological bloch bands in graphene superlattices. Proc. Natl Acad. Sci. USA 112, 10879–10883 (2015).

    ADS  MathSciNet  PubMed  PubMed Central  MATH  CAS  Google Scholar 

  47. 47.

    Zhang, Y.-H., Mao, D. & Senthil, T. Twisted bilayer graphene aligned with hexagonal boron nitride: anomalous hall effect and a lattice model. Phys. Rev. Res. 1, 033126 (2019).

    CAS  Google Scholar 

  48. 48.

    Xie, M. & MacDonald, A. H. Nature of the correlated insulator states in twisted bilayer graphene. Phys. Rev. Lett. 124, 097601 (2020).

    ADS  PubMed  CAS  Google Scholar 

  49. 49.

    Liu, S., Khalaf, E., Lee, J. Y. & Vishwanath, A. Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality. Phys. Rev. Res. 3, 013033 (2021).

    CAS  Google Scholar 

  50. 50.

    Liu, J. & Dai, X. Theories for the correlated insulating states and quantum anomalous hall effect phenomena in twisted bilayer graphene. Phys. Rev. B 103, 035427 (2021).

    ADS  CAS  Google Scholar 

  51. 51.

    Bultinck, N., Chatterjee, S. & Zaletel, M. P. Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene. Phys. Rev. Lett. 124, 166601 (2020).

    ADS  MathSciNet  PubMed  CAS  Google Scholar 

  52. 52.

    Scuri, G. et al. Electrically tunable valley dynamics in twisted wse 2/wse 2 bilayers. Phys. Rev. Lett. 124, 217403 (2020).

    ADS  PubMed  CAS  Google Scholar 

  53. 53.

    Schrade, C. & Fu, L. Spin-valley density wave in moiré materials. Phys. Rev. B 100, 035413 (2019).

    ADS  CAS  Google Scholar 

  54. 54.

    Zheng, Z. et al. Unconventional ferroelectricity in moiré heterostructures. Nature 588, 71–76 (2020).

    ADS  PubMed  CAS  Google Scholar 

  55. 55.

    Gu, J. et al. Dipolar excitonic insulator in a moire lattice. Preprint at (2021).

  56. 56.

    Zhang, Z. et al. Correlated interlayer exciton insulator in double layers of monolayer wse2 and moiré ws2/wse2. Preprint at (2021).

  57. 57.

    Jin, C. et al. Imaging of pure spin-valley diffusion current in ws2-wse2 heterostructures. Science 360, 893–896 (2018).

    ADS  PubMed  CAS  Google Scholar 

  58. 58.

    Li, T. et al. Quantum anomalous hall effect from intertwined moiré bands. Preprint at (2021).

  59. 59.

    Zhsng, Y., Devakul, T. & Fu, L. Spin-textured chern bands in ab-stacked transition metal dichalcogenide bilayers. Proc. Natl Acad. Sci. USA 118, e2112673118 (2021).

    Google Scholar 

  60. 60.

    Devakul, T., Crepel, V., Zhang, Y. & Fu, L. Dataset: Magic in transition metal dichalcogenide bilayers. Zenodo (2021).

Download references


We thank Kin Fai Mak, Jie Shan, Tingxin Li, and Shengwei Jiang for ongoing collaborations on MoTe2/WSe2, Bi Zhen and Constantin Schrade for previous collaborations on related topics. We thank Pablo Jarillo-Herrero, Kenji Yasuda, Cory Dean, Abhay Pasupathy, Qianhui Shi, Augusto Ghiotto, and En-Min Shih for helpful discussions. This work is primarily supported by the DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-SC0020149 (band structure calculation), DE-SC0018945 (theoretical modeling), and Simons Investigator award from the Simons Foundation (numerical analysis). L.F. is partly supported by the David and Lucile Packard Foundation.

Author information




T.D., V.C., Y.Z. and L.F. performed research, analyzed data, and wrote the manuscript.

Corresponding authors

Correspondence to Trithep Devakul or Liang Fu.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Devakul, T., Crépel, V., Zhang, Y. et al. Magic in twisted transition metal dichalcogenide bilayers. Nat Commun 12, 6730 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:

Further reading


By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.


Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing