Magic in twisted transition metal dichalcogenide bilayers

1In 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, the Tomonaga-Luttinger model for one-dimensional electron liquid,, and the Kitaev model for non-Abelian quantum spin liquid, 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 (TBG), experiments have discovered a variety of correlated electron states– facilitated by flat moiré bands.

2More 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–. 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”. The atomic potential depth and interatomic bonding are highly tunable by the choice of TMD materials, the twist angle and the displacement field,. Thus, TMD-based moiré materials provide a favorable platform for simulating idealized models in two dimensions.

3In 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.

4At 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. 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,.

5Due to spin-valley locking, 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.

6In order to obtain accurate moiré band structures, we perform large-scale density functional theory calculations with the SCAN+rVV10 van der Waals density functional, 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 studies and consistent with recent works,,.

7At 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., 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 byand as its time-reversal conjugate.

8The 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,, 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 vector. In the lowest harmonic approximation, ΔT(r) and Δ1,2(r) are sinusoids that interpolate between MM, MX and XM regions:where gj are (j − 1)π/3 counter-clockwise rotations of the moiré reciprocal lattice vector, and a0 is the monolayer lattice constant. Up to an overall energy scale, the continuum model depends only on the dimensionless parameters, w/V and ψ.

9From 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. 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).

10We remark that different approaches,, 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). 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.

11Using the new continuum model parameters established above, along with the lattice constant a0 = 3.317 Å and the effective mass m* = 0.43me,, 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,, as well as the (direct or indirect) band gaps εij between pairs of bands (i, j) = (1, 2) and (2, 3),, 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∘.

12First, for θ < θ1, the top two bands are well separated from the rest of the spectrum, and carry opposite Chern number. 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. 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.

13At θ = θ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 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.

14For θ < θ2, especially near θ2 where the first band is more dispersive, the spin Chern number 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, giving rise to a double quantum spin Hall state with two sets of counter-propagating spin-polarized edge modes at n = 4.

15We 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.

16Analytic 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, the effective mass mγ should diverge near the crossover. Let to be the angle at which. 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, and n0 is the integer (mod 6) which maximizes (n0 = 1 for WSe2 parameters). We find that Eq. (4) provides a decent estimate for θm in the cases considered. In WSe2, we have, compared to θm = 1.43∘ at which the bandwidth is minimized.

17Next, we turn to the Berry curvature of the top band, shown in Fig. 2b. In all cases, the first band has Chern number; however, the distribution changes very drastically as θ is varied. At θ ≳ 2∘, is peaked around the band minimum at γ. At θ ≲ 1∘, is instead sharply peaked at the κ± points. Near the crossover region, the distribution of shifts from γ to κ±, and can become very evenly distributed. We find that is most evenly distributed near θ = 1.67∘, shown in Fig. 2b, where becomes almost uniform in the Brillouin zone. The uniform distribution of is reminiscent to that of Landau levels. Time-reversal symmetry forces the corresponding spin-↓ bands from the −K valley to have opposite Chern number.

18We 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.

19Recall 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 of the continuum Hamiltonian (1) for each k in the mBZ, we first construct a superposition of the first two (n = 1, 2) eigenstates, using a 2 × 2 unitary matrix, which maximizes the layer polarization at every k:where is the projector on to the top/bottom layer, so that is chosen to mostly consist of states in the top layer, and similarly for on the bottom layer. This uniquely specifies 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. 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: is mostly in the top layer, and vice versa for.

20It 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 regions.

21The tight binding Hamiltonian is found to bewhere 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 ≡ ∣t2∣eiϕ 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.

22For 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. 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 strength,.

23In 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 with σ = ± are exact eigenstates of our interacting model, because the spin-orbit coupling in 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 exactly corresponds to Haldane’s model for QAH insulator.

24To 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, 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, 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∘.

25To 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. 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.

26For 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.

27We 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.

28It 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,, where the alignment with hBN substrate is likely the origin of valley Chern number and both spin and valley degeneracy are lifted due to repulsive interaction in the flat band–. 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.

29Another great advantage of twisted TMD bilayer is their high degree of tunability, in particular with respect to applied electric fields,,. Due to the layer polarization of the Wannier basis states, the displacement field can be modeled as a sublattice symmetry breaking term, 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 = 10∣t1∣, 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 unit cell in order to describe all ordered phases of the model.

30At 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 (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). 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).

31Beyond 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 t2 ≲ t1 ≪ Δ, U relevant for our system, their coupling is described by an effective XXZ model with Dzyaloshinskii−Moriya (DM) interactionswhich 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 bywith. In Eq. 8, we have separated exchange terms coming from different physical processes. The first ones arise from nearest neighbor tunneling on the triangular A sublattice, while the others originate from loop-exchange on the honeycomb lattice that do not involve any double occupancy.

32For twist angles θ ≲ 1∘, t2 ≪ t1 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,,. 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 density wave—has been proposed for the correlated insulating state at n = 1 in twisted bilayer WSe2 at θ ~ 4∘–5∘.

33For 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 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. 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.

34Finally, we comment on the effect of nearest neighbor repulsion V∑〈i, j〉ninj 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 ferroelectricity,. The long-range component of interactions can be controlled by screening from nearby metallic layers. Multiple recent experiments, 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.

35Our 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 transport, which can be detected by optical spin injection and spatial-temporal mapping recently developed for TMD bilayers. 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.

36In 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 field,. 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.

37Looking 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.