## Abstract

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 WSe_{2}, 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.

## Introduction

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 system^{1}, the Tomonaga-Luttinger model for one-dimensional electron liquid^{2,3}, and the Kitaev model for non-Abelian quantum spin liquid^{4}, 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 graphene^{5} (TBG), experiments have discovered a variety of correlated electron states^{6,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 lattice^{13,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 angle^{13} and the displacement field^{20,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 model^{33}.

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 insulator^{34}. 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 studies^{23,35}.

Due to spin-valley locking^{36}, monolayer TMDs such as WSe_{2} and MoTe_{2} 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.

## Results

In order to obtain accurate moiré band structures, we perform large-scale density functional theory calculations with the SCAN+rVV10 van der Waals density functional^{37}, which captures the intermediate-range vdW interaction through its semilocal exchange term. Focusing on twisted bilayer WSe_{2}, 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 WSe_{2} are found to come from the ±*K* valley (shown in Fig. 1c), as opposed to the Γ valley in previous computational studies^{38} and consistent with recent works^{35,39,40}.

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

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 **d**_{0}. 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 **a**_{1,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 vector^{20}. In the lowest harmonic approximation, Δ_{T}(**r**) and Δ_{1,2}(**r**) are sinusoids that interpolate between MM, MX and XM regions^{14}:

where **g**_{j} 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 *a*_{0} 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 reported^{14}. 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 approaches^{31,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 *C*_{3z} eigenvalue at the high symmetry momenta *γ*, *κ*_{±}, which determine their Chern number (mod 3)^{41}. The *C*_{3z} 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\).

Using the new continuum model parameters established above, along with the lattice constant *a*_{0} = 3.317 Å^{42} and the effective mass *m*^{*} = 0.43*m*_{e}^{43,44}, we calculate the band structure of twisted bilayer WSe_{2}, *E*_{i}(**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^{∘}.

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^{∘}, *E*_{1}(**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 *E*_{1}(**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

*E*

_{1}(

**) 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)**

*γ*where \({{{{{{{{\mathcal{E}}}}}}}}}_{n}=2w\cos (\pi n/3)+2V\cos (2\pi n/3-\psi )\), and *n*_{0} is the integer (mod 6) which maximizes \({{{{{{{{\mathcal{E}}}}}}}}}_{n}\) (*n*_{0} = 1 for WSe_{2} parameters). We find that Eq. (4) provides a decent estimate for *θ*_{m} in the cases considered. In WSe_{2}, 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 singularity^{23}.

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**:

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 *a*_{W} ≈ 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 *n*th nearest neighbor hopping matrix elements *t*_{n} 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 regions^{14}.

The tight binding Hamiltonian is found to be

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 *t*_{1} is real, while *t*_{2} ≡ ∣*t*_{2}∣*e*^{iϕ} is complex, and ⋯ contain longer range hopping terms. We find that ∣*t*_{n}∣ 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 TB_{n} including up to *t*_{n} hopping terms, compared to that of the continuum model. For *θ* ≲ 1^{∘}, TB_{2} 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* ~ *e*^{2}/(*ϵ**a*_{W}) ≈ 70 meV at *θ*_{m}, using a realistic relative dielectric constant *ϵ* = 4 and *a*_{W} = 5nm, which is significantly larger than the tight binding parameters *t*_{n}. Therefore at such small twist angle, twisted WSe_{2} homobilayers are in the strong-coupling regime, in contrast with *θ* ~ 4^{∘}–5^{∘} where the bandwidth is comparable to the interaction strength^{23,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^{∘}.

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 subspace^{45}. 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 WSe_{2}, 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 TBG^{11,12}, where the alignment with hBN substrate is likely the origin of valley Chern number^{46} and both spin and valley degeneracy are lifted due to repulsive interaction in the flat band^{47,48,49,50,51}. Due to the presence of *S**U*(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 fields^{19,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 *s*_{i} 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* = 10∣*t*_{1}∣, 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 Δ ≲ 2*t*_{1} (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 *t*_{n≥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 (FM_{z})^{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^{∘} AFM_{xy}), or ferromagnetically align in the *x**y* plane (FM_{xy}). 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 *t*_{2} ≲ *t*_{1} ≪ Δ, *U* relevant for our system, their coupling is described by an effective XXZ model with Dzyaloshinskii−Moriya (DM) interactions

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

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^{∘}, *t*_{2} ≪ *t*_{1} 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 TMDs^{24,26,53}. It was shown to yield an antiferromagnetic phase that the small residual DM interaction pins in the *x**y* plane. This is the origin of the AFM_{xy} phase observed in Fig. 5b. We also note that the weak-coupling version of AFM_{xy} phase—an intervalley-coherent \(\sqrt{3}\times \sqrt{3}\) density wave^{23}—has been proposed for the correlated insulating state at *n* = 1 in twisted bilayer WSe_{2} at *θ* ~ 4^{∘}–5^{∘}^{35}.

For larger twist angles, *t*_{2} 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 FM_{xy} phase is favored as shown in Fig. 5b. The competition between AFM_{xy} and FM_{xy} 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 Δ = 5*t*_{1}, 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 *V*∑_{〈i, j〉}*n*_{i}*n*_{j} to our HF phase diagram. This term favors the layer polarized phases, such as the FM_{xy} and AFM_{xy} 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 ferroelectricity^{20,54}. The long-range component of interactions can be controlled by screening from nearby metallic layers. Multiple recent experiments^{55,56} on twisted WSe_{2} homobilayers in the presence of a nearby WSe_{2} 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.

## Discussion

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: FM_{z}, FM_{xy}, and AFM_{xy}. Despite being electrically insulating, the *x**y*-ordered Mott insulators support coherent magnon transport^{23}, which can be detected by optical spin injection and spatial-temporal mapping recently developed for TMD bilayers^{57}. 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 MoTe_{2}/WSe_{2} at the filling of *n* = 1 tuned by displacement field^{58,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 database^{60}. Additional data related to this paper are available from the corresponding authors upon reasonable request.

## References

- 1.
Hubbard, J. Electron correlations in narrow energy bands.

*Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci.***276**, 238–257 (1933). - 2.
Tomonaga, S.-I. Remarks on Bloch’s method of sound waves applied to many-Fermion problems.

*Prog. Theor. Phys.***5**, 544–569 (1950). - 3.
Luttinger, J. M. An exactly soluble model of a many fermion system.

*J. Math. Phys.***4**, 1154–1162 (1963). - 4.
Kitaev, A. Anyons in an exactly solved model and beyond.

*Ann. Phys.***321**, 2–111 (2006). - 5.
Bistritzer, R. & MacDonald, A. H. Moiré bands in twisted double-layer graphene.

*Proc. Natl Acad. Sci. USA***108**, 12233–12237 (2011). - 6.
Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices.

*Nature***556**, 80–84 (2018). - 7.
Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices.

*Nature***556**, 43–50 (2018). - 8.
Yankowitz, M. et al. Tuning superconductivity in twisted bilayer graphene.

*Science***363**, 1059–1064 (2019). - 9.
Lu, X. et al. Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene.

*Nature***574**, 653–657 (2019). - 10.
Cao, Y. et al. Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene.

*Nature***583**, 215–220 (2020). - 11.
Sharpe, A. L. et al. Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene.

*Science***365**, 605–608 (2019). - 12.
Serlin, M. et al. Intrinsic quantized anomalous hall effect in a moiré heterostructure.

*Science***367**, 900–903 (2020). - 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). - 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). - 15.
Tang, Y. et al. Simulation of hubbard model physics in wse2/ws2 moiré superlattices.

*Nature***579**, 353–358 (2020). - 16.
Regan, E. C. et al. Mott and generalized wigner crystal states in wse2/ws2 moiré superlattices.

*Nature***579**, 359–363 (2020). - 17.
Shabani, S. et al. Deep moiré potentials in twisted transition metal dichalcogenide bilayers.

*Nat. Phys.***17**, 720–725 (2021). - 18.
Jin, C. et al. Stripe phases in wse2/ws2 moiré superlattices.

*Nat. Mater.***20**, 940–944 (2021). - 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). - 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). - 21.
Slagle, K. & Fu, L. Charge transfer excitations, pair density waves, and superconductivity in moiré materials.

*Phys. Rev. B***102**, 235423 (2020). - 22.
Xu, Y. et al. Correlated insulating states at fractional fillings of moiré superlattices.

*Nature***587**, 214–218 (2020). - 23.
Bi, Z. & Fu, L. Excitonic density wave and spin-valley superfluid in bilayer transition metal dichalcogenide.

*Nat. Commun.***12**, 642 (2021). - 24.
Pan, H., Wu, F. & Das Sarma, S. Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer wse

_{2}.*Phys. Rev. Res.***2**, 033087 (2020). - 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.
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.
Padhi, B., Chitra, R. & Phillips, P. W. Generalized wigner crystallization in moiré materials.

*Phys. Rev. B***103**, 125146 (2021). - 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). - 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). - 30.
Zhang, Z. et al. Flat bands in twisted bilayer transition metal dichalcogenides.

*Nat. Phys.***16**, 1093–1096 (2020). - 31.
Magorrian, S. J. et al. Multifaceted moiré superlattice physics in twisted wse

_{2}bilayers.*Phys. Rev. B***104**, 125440 (2021). - 32.
Tang, H., Carr, S. & Kaxiras, E. Geometric origins of topological insulation in twisted layered semiconductors.

*Phys. Rev. B***104**, 155415 (2021). - 33.
Kane, C. L. & Mele, E. J. Quantum spin hall effect in graphene.

*Phys. Rev. Lett.***95**, 226801 (2005). - 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). - 35.
Wang, L. et al. Correlated electronic phases in twisted bilayer transition metal dichalcogenides.

*Nat. Mater.***19**, 861–866 (2020). - 36.
Xiao, D., Liu, G.-B., Feng, W., Xu, X. & Yao, W. Coupled spin and valley physics in monolayers of mos

_{2}and other group-vi dichalcogenides.*Phys. Rev. Lett.***108**, 196802 (2012). - 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). - 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). - 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.
Kundu, S., Naik, M. H., Krishnamurthy, H.R. & Jain, M. Flat bands in twisted bilayer wse _2 with strong spin−orbit interaction. Preprint at https://arxiv.org/abs/2103.07447 (2021).

- 41.
Fang, C., Gilbert, M. J. & Bernevig, B. A. Bulk topological invariants in noninteracting point group symmetric insulators.

*Phys. Rev. B***86**, 115112 (2012). - 42.
Mounet, N. et al. Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds.

*Nat. Nanotechnol.***13**, 246–252 (2018). - 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). - 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). - 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). - 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). - 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). - 48.
Xie, M. & MacDonald, A. H. Nature of the correlated insulator states in twisted bilayer graphene.

*Phys. Rev. Lett.***124**, 097601 (2020). - 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). - 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). - 51.
Bultinck, N., Chatterjee, S. & Zaletel, M. P. Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene.

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

*Phys. Rev. Lett.***124**, 217403 (2020). - 53.
Schrade, C. & Fu, L. Spin-valley density wave in moiré materials.

*Phys. Rev. B***100**, 035413 (2019). - 54.
Zheng, Z. et al. Unconventional ferroelectricity in moiré heterostructures.

*Nature***588**, 71–76 (2020). - 55.
Gu, J. et al. Dipolar excitonic insulator in a moire lattice. Preprint at https://arxiv.org/abs/2108.06588 (2021).

- 56.
Zhang, Z. et al. Correlated interlayer exciton insulator in double layers of monolayer wse2 and moiré ws2/wse2. Preprint at https://arxiv.org/abs/2108.07131 (2021).

- 57.
Jin, C. et al. Imaging of pure spin-valley diffusion current in ws2-wse2 heterostructures.

*Science***360**, 893–896 (2018). - 58.
Li, T. et al. Quantum anomalous hall effect from intertwined moiré bands. Preprint at https://arxiv.org/abs/2107.01796 (2021).

- 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). - 60.
Devakul, T., Crepel, V., Zhang, Y. & Fu, L. Dataset: Magic in transition metal dichalcogenide bilayers. Zenodo https://doi.org/10.5281/zenodo.5607764 (2021).

## Acknowledgements

We thank Kin Fai Mak, Jie Shan, Tingxin Li, and Shengwei Jiang for ongoing collaborations on MoTe_{2}/WSe_{2}, 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

### Affiliations

### Contributions

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

### Corresponding authors

## 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 http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Devakul, T., Crépel, V., Zhang, Y. *et al.* Magic in twisted transition metal dichalcogenide bilayers.
*Nat Commun* **12, **6730 (2021). https://doi.org/10.1038/s41467-021-27042-9

Received:

Accepted:

Published:

DOI: https://doi.org/10.1038/s41467-021-27042-9

## Further reading

## Comments

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.