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.

Density functional theory calculation and continuum model -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 (AB stacking can be viewed as a 180 • rotation of top layer). Within a local region of a twisted bilayer, the atom configuration is identical to that of an untwisted bilayer, where one layer is laterally shifted relative to the other layer by a corresponding displacement vector d 0 . For this reason, the moiré band structures of twisted TMD bilayers can be constructed from a family of untwisted bilayers at various d 0 , all having 1 × 1 unit cell. Our analysis thus starts from untwisted bilayers.
In particular, d 0 = 0, (−a 1 + a 2 ) /3, (a 1 + a 2 ) /3, where a 1,2 is the primitive lattice vector for untwisted bilayers, correspond to three high-symmetry stacking configurations of untwisted TMD bilayers, which we refer to as MM, XM, MX. In MM (MX) stacking, the M atom on the top layer is locally aligned with the M (X) atom on the bottom layer, likewise for XM. The bilayer structure in these stacking configurations is invariant under three-fold rotation around the z axis.
Density functional calculations are performed using generalized gradient approximation with SCAN+rVV10 Van der Waals density functional, as implemented in the Vienna Ab initio Simulation Package. Pseudopotentials are used to describe the electron-ion interactions. We first construct the zero-twisted angle WSe2/WSe2 bilayer at MM and MX lateral configurations with vacuum spacing larger than 20Å to avoid artificial interaction between the periodic images along the z direction. The structure relaxation is performed with force on each atom less than 0.01 eV/A. We use 12 × 12 × 1 for structure relaxation and self-consistent calculation. The more accurate SCAN+rVV10 van der Waals density functional gives the relaxed layer distances as 6.6Å and 7.15Å for MX and MM stacking structures, respectively. By calculating the work function from electrostatic energy of converged charge density, we plot in Fig. 1 the band structure of MM and MX-stacked bilayers, with reference energy E = 0 chosen to be the absolute vacuum level.
To compare with moiré band structure, we further construct the commensurate structures with twist angle 5.08 • with 762 atoms, and perform large scale DFT cal- culation to relax the moiré superlattice. The spacial dependent layer distance profile has been discussed in the main text, and fits nicely with the untwisted MM and MX stacking structures. The moiré band structure with spin orbital coupling agrees remarkably well with the continuum model derived from untwisted calculations, and the effective mass is chosen as 0.38m e , close to the experimental value 0.4∼0.45m e . To reveal the real space localization of moiré flat band, the γ point Kohn-Sham wavefunction for the top-most moiré valence band is calculated and projected to bottom and top layer, respectively. As shown in Fig. 2, we find that the charge are localized at MX and XM region, consistent with previous Wannier function. For MoTe 2 , the Γ pockets appears to be the valence band maximum at MX stacking structure when using lattice constant 3.52Å from bulk structures (The relaxed layer distance is 6.95Å, close to the bulk layer distance 6.98Å in 2H structures). Experimentally, AB stacked bilayer MoTe 2 has been shown to be a indirect band gap semiconductor from optical measurements [1,2]. We thus believe Γ pockets would enter in the moiré valence bands for twisted homobilayer MoTe 2 even for AA stacking. At θ = 5.08 • , we calculate the moiré band structure of fully relaxed AA stacked homobilayer MoTe 2 . As shown in Fig. 3, the narrow Dirac like Γ valley moire bands intersect with topmost dispersive K valley moire bands. In transport measurements under zero displacement field, these two type of charge carriers would both enter, causing further complications.

B. Supplementary Note 2
Moiré topology criteria -We determine the Chern number of top two moiré valence bands from their C 3z eigenvalues. Computationally, the Kohn-Sham wavefunctions in plane-wave basis are extracted at three C 3z invariant momentum γ, κ, and κ with energy cutoff as E cut = 100 eV and 2 (k + G) 2 /2m e < E cut . And we compute their C 3z eigenvalue as: where σ, σ =↑ or ↓, Ψ σ nk |C 3z |Ψ σ nk calculated from the transformation relation of plane-waves, and S σσ (C 3z ) is the C 3z spin rotation matrix.
As there is a valley time reversal symmetry within the moire superlattice, the energy bands between γ, κ and κ are enforced to be double degenerate from K and −K valley. Therefore, to determine the valley index of moiré bands, we add a small gating field D = 0.005 V/nm to introduce a less than 2 meV splitting at κ and κ . The K and −K indexes are then diagnosed by spin resolved wavefunction. At γ point where bands are still double degenerate with gating field, we extracted the C 3z eigenvalues and valley indexes from finite momentum perturbation along γ to κ direction.

C. Supplementary Note 3
Continuum model at higher angles - Figure 4 show the band structure for the continuum model at larger twist angles than discussed in the main text, as well as the band gaps and bandwidth. The bandwidth of the first band W 1 increases rapidly for θ > θ 1 ≈ 1.5 • , reaching W ∼ 100meV for θ ∼ 5 • -6 • . The bandwidth of the second band W 2 has a minima at θ ≈ 2.2 • , and remains quite small throughout the angle range Analytic magic angle -In this section, we analytically estimate the magic angle by consider the effective mass of the top band near γ.
At k = γ, when w = V = 0, the highest energy states are sixfold degenerate plane wave states at energy e 0 = − |γ−κ±| 2 , which is separated from the next set of eigenstates at . Incorporating the effect of w and V will introduce coupling between these states at e 0 , as well as with those at e 1 and beyond. For the purposes obtaining an analytic estimate for the magic angle, we include only this set of 6 highest energy states. This approximation is valid for w, V e 0 −e 1 : as a rough estimate of the validity, e 0 − e 1 ≈ 130meV at θ = 1.5 • is an order of magnitude greater than (V, w) = (9, 18)meV for WSe 2 .
These six states consist of the plane waves at k j for j = 1 . . . 6, where and R θ is a θ counter-clockwise rotation matrix. The continuum Hamiltonian acts on these states as where |j ± 6 ≡ |j , and · · · contain states outside of {|j }, which we have ignored due to the large energy gap e 1 −e 0 . Utilizing the emergent C 6 symmetry which sends |j → |j + 1 , H ↑ is easily diagonalized (in the restricted subspace) by eigenstates with eigenvalue e 0 + E n , where E n = 2w cos(πn/3) + 2V cos(2πn/3 − ψ).
In the absence of any accidental degeneracy, there is a unique highest energy eigenstate |n 0 . Next, in order to estimate when γ transitions from band maxima to minima, we consider a small momentum shift about γ. The C 3 symmetry of the microscopic model restricts the lowest order momentum dependence to be simply 1 2mγ |k − γ| 2 for some effective mass m γ . Thus, we consider a momentum shift x away from γ and work perturbatively in . Expanding the energy E = E (0) + 2 E (2) , we have that E (2) = 1/(2m γ ), and the magic angleθ m is where E (2) changes sign.
We expand H ↑ as H ↑ = H 0 + H 1 + 2 H 2 , which act on |j as and H 2 = − 1 2m * . We havê hence, acting on |n , Applying standard perturbation theory in the |n basis, the second order correction to the energy is Solving for where E (2) = 0 gives the magic angle condition from the main text, Localized Wannier functions and tight binding model -In this section, we provide additional details on the construction of the Wannier states shown in the main text. We numerically diagonalize the Hamiltonian H ↑ on an N k = 30×30 discretization of the moiré Brillouin zone, corresponding to a system of 30×30 moiré unit cells. For each k, the eigenstates {|φ n,k } n is obtained in the plane wave basis, expressed as |φ nk = g,l c (nk) gl |k + g, l , where g are reciprocal moiré lattice vectors, l = ± is the layer, and c (nk) gl are complex numbers obtained from exact diagonalization. The diagonalization is done by truncating to some maximum |k + g| which is sufficiently far from the valence band top at κ ± so as to have negligible effect on the first few moiré bands.
Constructing Wannier states for the top two bands involves finding the appropriate 2×2 unitary matrices U (k) nm at each k. For convenience, we shall suppress all k labels. We parameterize the unitary matrix as U = e iψ1 cos(χ/2) e iψ1+iδ sin(χ/2) e iψ2−iδ sin(χ/2) −e iψ2 cos(χ/2) , where χ ∈ [0, π], and δ, ψ 1 , ψ 2 ∈ [0, 2π]. As stated in the main text, U is chosen to maximize the layer polarization of |φ n = m U nm |φ m , where P ± projects to the ± layer. Using P − = 1 − P + , this can be simplified to where P nm = φ n |P + |φ m = P * mn . It is straightforward to maximize P with respect to χ and δ. First, setting ∂ δ P = 0 gives δ = arg P 21 + sπ, where s ∈ {0, 1}. Then, setting ∂ χ P = 0 gives (notice the range of arctan is taken to be [0, π]). The two solutions for s = 0, 1 correspond to the maximum and minimum of P . To find the maximum, we simply compare the value of P for these two solutions and take the larger one.
To determine the phase factor ψ 1,2 , we take the real space component of the wavefunction at the MX/XM stacking regions. Define r 1 = (a 0 /2θ)[−1/ √ 3, 1] and r 2 = (a 0 /2θ)[1/ √ 3, 1] to be the locations of the XM and MX stacking region. We then take the phase factors ψ 1,2 such that r n , l n |φ 1 is real and positive, where l 1(2) = +(−) (here, |r, l is the usual real space state at r on layer l, satisfying r, l|k + g, l = e i(k+g)·r δ ll ). This fixes the phase factors ψ n = − arg r n , l n |φ n .
Finally, the Wannier functions localized at R + r n , for a moiré lattice vector R, are defined The Wannier centers R + r n form a moiré honeycomb lattice on the XM and MX stacking regions. The hopping matrix elements between Wannier states on the site R + r n and R + r n can be obtained from the continuum Hamiltonian and U matrices. We have where E m (k) is the energy of |φ mk under H ↑ . Using this method, we compute the hopping matrix elements between nth nearest neighbor sites up to n = 5 as illustrated in Fig 5, and shown in the main text. We find that t n are mostly real with the exception of t 2 , which we define (see Fig 5) to be the 2nd nearest neighbor hopping term along right-turning path. From the symmetries of the model, it follows that the hopping term along a leftturning path is t * 2 . The −K valley produces a similar tight binding model of spin-↓ electrons related to this one by time reversal symmetry. As a final remark, we note that the Wannier states obtained through this method are not maximally localized, thus our tight binding model has the potential to be further improved.

F. Supplementary Note 6
Effective Spin Hamiltonian -The derivation of an effective spin Hamiltonian from a perturbative expansion of Eq. 6 in the limit |t 1 |, |t 2 | ∆, U is rather complicated when no clear scale separation exists between ∆ and U . Fortunately in our model, the interaction scale U largely dominates over the tunneling terms of H TB , and all experimentally relevant applied displacement fields ∆. As a consequence, we may first eliminate all doubly occupied configurations, and then use t 1 ∆ to obtain the effective spin Hamiltonian given in the main text.

Strong coupling expansion
Let us first fix some notations. We write the interacting Kane-Mele model as where T u gathers all single-particle terms which change the number of doubly occupied by u. To get their explicit expression, we gather in T δ,u the nearest neighbor tunneling processes changing the number of occupied A site by δ and the number of doubly occupied sites by u, where δ = ± and u = 0, ±. A similar definition T (2) 0,u applies for next-nearest neighbor hoppings. These definition yield Elimination of doubly occupied configurations is achieved by standard perturbation theory techniques, see for instance Ref. [3], and gives that we have recast in a form more amenable to perturbative expansion in ∆ by introducing 0,+ + T 0,+ , Similar to the previous case,T δ changes the number of occupied A sites by δ. When their potential energy is large t 1 ∆, these A-sites can be eliminated in a similar way as before, see Ref. [4] for details. This leads to the following effective Hamiltonian in the physical subspace with no double occupations and no occupied A sites: This expression can be greatly simplified by noticing that the direct action of T ±,± and T 0,0 is zero on the physical space with no double occupancy and no occupied A sites. Furthermore, all terms of order U −2 should be discarded for consistency with our first expansion Eq. 19. Using these insights, we obtain where we have isolated the contribution 0,+ + hc Supplementary Figure 6: Multi-tunneling processes appearing in the perturbation expansion for large displacement field and large interaction scale. a) Effective tunnelingt = t 2 + t 2 1 /∆ on the B triangular lattice after elimination of the A site with large potential energy. b) Spin exchange ∝t 2 /U between localized spins. c) Additional loop exchange ∝ t 2 1 t 2 /U process which does not involve any double occupations. The choice of t 2 /t * 2 in the diagrams applies to spin-up electrons. Grey circles denote doubly occupied B sites, while arrows represents the spin of the holes on A sites. that we from now on neglect since its magnitude |t 1 t 2 | 2 /(U ∆ 2 ) is much smaller than the other terms appearing in H U ∆ tot when |t 2 | |t 1 |, as is the case for our effective tightbinding model (see I E).
Discarding an overall energy shift, we can rewrite the Hamiltonian derived above as The first line simply describes the virtual occupation of a doubly excited state using the bare or an effective nextnearest tunneling, while the second line corresponds to a particle exchange process where particle turn around a small loop of the honeycomb lattice to avoid double occupations. They are sketched in Fig. 6. It is worth noting here that the terms obtained in H S are simply those one would have guessed without decoupling the U and ∆ expansions. This is however only due to the facts that the operatorsT ±2 are zero in the physical space considered in our model. This is why we have preferred the more lengthy, but more reliable, method presented here.

Rewriting with spin operators
Introducing the effective tunnelingt = t 2 + t 2 1 /∆, we can rewrite the first line of Eq. 24 as which corresponds to the processes sketched in Fig. 6a and 6b. With a similar calculation, we can recast the second line as which is sketched in Fig. 6c. Putting everything together, we obtain an XXZ model with Dzyaloshinskii-Moriya interactions for the spin localized on the triangular B lattice: where the coefficients read This is the form used in the main text.

Classical phase diagram
The phases appearing for large displacement fields in our Hartree-Fock study (Fig. 5 in the main text) can be rationalized by studying the classical version of Eq. 27. Treating the spins as classical vectors of fixed lengths and imposing the 'weak constraint' of Luttinger-Tisza [5], the spin model can be exactly solved in momentum space.
We find spin eigenstates with spin polarized along z with the dispersion E z (k) = J C(k), and two other with spins pined in the xy plane where we have introduced the two auxiliary functions cos(k · a j ), S(k) = 3 j=1 sin(k · a j ), (31) with a 3 = −(a 1 + a 2 ). The competition between the out of plane and in-plane modes was studied beyond the classical approximation with the ED spin-wave analysis (see main text). Thus, we discard the E z eigenstates and look for spin patterns that minimize E ± xy . Depending on J ⊥ and D, the minimum of E − xy is either at the γ point or at κ ± , which respectively describe the FM xy and AFM xy phases of the main text, after reintroduction of Luttinger-Tisza's 'strong constraint' [5]. Comparison of the energies at these three momenta give the conditions FM xy : |D| + √ 3J ⊥ < 0, AFM xy : |D| + √ 3J ⊥ > 0.