Interaction-driven topological phase diagram of twisted bilayer MoTe 2

Twisted bilayer MoTe 2 is a promising platform to investigate the interplay between band topology and many-body interaction. We present a theoretical study of its interaction-driven quantum phase diagrams based on a three-orbital model, which can be viewed as a generalization of the Kane-Mele-Hubbard model with one additional orbital and long-range Coulomb repulsion. We predict a cascade of phase transitions tuned by the twist angle θ . At the hole filling factor ν = 1 (one hole per moir´e unit cell), the ground state can be in the multiferroic phase with coexisting spontaneous layer polarization and magnetism, the quantum anomalous Hall phase, and finally the topologically trivial magnetic phases, as θ increases from 1 . 5 ◦ to 5 ◦ . At ν = 2, the ground state can have a second-order phase transition between an antiferromagnetic phase and the quantum spin Hall phase as θ passes through a critical value. The dependence of the phase boundaries on model parameters such as the gate-to-sample distance, the dielectric constant, and the moir´e potential amplitude is examined. The predicted phase diagrams can guide the search for topological phases in twisted transition metal dichalcogenide homobilayers.

In this paper, we present a theoretical study of the interaction-driven quantum phase diagram of tMoTe 2 , where the two layers are rotated by an angle θ relative to the rhombohedral-stacked bilayer.As predicted in Ref. 2, low-energy moiré valence bands of tMoTe 2 originate from arXiv:2305.01006v3 [cond-mat.mes-hall]15 Nov 2023 ±K valleys of monolayer MoTe 2 , where ±K refer to the monolayer Brillouin zone corners.The intralayer moiré potentials confine low-energy states to two different locations [i.e., A and B sites of the moiré superlattice shown in Fig. 1(a)] in the bottom and top layers, where the two locations exactly correspond to the two sublattices of a honeycomb lattice.The rotation between the two layers generates layer-dependent momentum shifts in the kinetic energy, which induces effective valley-contrasting and sublattice-dependent fluxes on the honeycomb lattice model.Because of the lattice geometry and the flux pattern, the first two moiré valence bands realize the Haldane model [45] within one valley and the Kane-Mele model [46,47] when both valleys are considered [2].Moreover, the moiré bands have narrow bandwidth, and electron interaction effects are enhanced.Therefore, tMoTe 2 can serve as a model system to study the interplay between band topology and many-body interaction.However, tMoTe 2 has not been studied experimentally until very recently [21][22][23].Two independent experiments reported signatures of integer and fractional quantum anomalous Hall states in tMoTe 2 with θ ≈ 3.7 • [22] and 3.4 • [23], respectively.This exciting experimental progress makes the theoretical prediction of Ref. 2 to be of experimental relevance, and calls for thorough studies of the many-body interaction effects in tMoTe 2 .
The mapping between the first two moiré bands and the Kane-Mele model allows to study the interaction effects within a lattice model, which should be a generalized Kane-Mele-Hubbard model since the realistic interaction is the long-range Coulomb repulsion.However, such a mapping is limited to the small θ regime of θ < θ * 1 , where θ * 1 is estimated to be about 1.74 1 , there is a topological band inversion between the second and the third bands.To describe the low-energy moiré bands in a unified manner regardless of the exact values of θ, we construct a three-orbital model based on the Wannier states of the first three bands.The threeorbital model faithfully captures the energy, symmetry, and topology of the first three bands.The interacting model is then obtained by projecting the Coulomb interaction onto the Bloch-like states derived from the Wannier orbitals.The interacting Hamiltonian is formulated in momentum space and implicitly includes all microscopic interactions such as density-density interaction and Hund's coupling.
We calculate the quantum phase diagrams of the interacting three-orbital model at hole filling factors ν = 1 and ν = 2, where ν counts the number of doped holes per moiré unit cell.The phase diagrams are obtained by comparing the energy of multiple competing states within self-consistent Hartree-Fock approximation.We construct the phase diagram as a function of θ and examine its dependence on the gate-to-sample distance d, the dielectric constant ϵ, and the moiré potential amplitude V .We predict a cascade of phase transitions tuned by θ in the range between 1.5 • and 5 • .At ν = 1, the ground state is in the multiferroic phase with spontaneous layer polarization and magnetism for small θ, the quantum anomalous Hall phase with a Chern number of 1 for intermediate θ, and topologically trivial magnetic phases for large θ.At ν = 2, the ground state is in interactiondriven magnetic phases for small θ and turns into the quantum spin Hall phase for large θ.We also discuss how the phase boundaries depend on parameters such as d, ϵ, and V .The predicted phase diagrams demonstrate the remarkable richness of tMoTe 2 , and are expected to guide further experimental and theoretical studies.
There are several benefits of doing the many-body calculation using the three-orbital model compared to the full continuum model.First, the constructed Wannier states provide a clear real-space picture of the low-energy degrees of freedom and valuable intuition for studying interaction physics.The localized Wannier orbitals allow us to devise various mean-field states with different layer polarization, magnetization, and topology.Second, the symmetry and topology of electronic states become more apparent in the Wannier basis.Third, the three-orbital model can be further studied using techniques (e.g., exact diagonalization) beyond mean-field theory.
The paper is organized as follows.In Sec.II, we construct the noninteracting as well as the interacting threeorbital model.Construction of the Wannier states is informed by symmetry analysis of the moiré bands, which is detailed in Appendix A. In Sec.III, we present the phase diagrams at ν = 1 and ν = 2.Additional discussion on the phase diagram and the model is given, respectively, in Appendices B and C. In Sec.IV, we conclude with a summary and a discussion.Our theoretical results are compared with the recent experimental findings, and future directions are discussed.

A. Moiré Hamiltonian
The moiré superlattices of tMoTe 2 have D 3 pointgroup symmetry generated by a threefold rotation C 3z around the out-of-plane ẑ axis and a twofold rotation C 2y around the in-plane ŷ axis that exchanges the bottom (b) and top (t) layers.As illustrated in Fig. 1(a), the moiré superlattices have three distinct high-symmetry locations labeled by O, A, and B, where the latter two positions are related by the C 2y symmetry.
The single-particle moiré Hamiltonian for valence band states in tMoTe 2 has been constructed in Ref. 2 and is given by, Here τ = ± represents ±K valleys and also spin up (↑) and down (↓), since spin and valley are locked for valence band states in TMD [7].The moiré Hamiltonian H τ in Eq. ( 1) is represented by a 2 × 2 matrix in the layerpseudospin space.The diagonal terms in H τ describe an electron moving in the layer-dependent potential ∆ ± (r) with an amplitude V and phase parameters ±ψ, and the off-diagonal terms account for the interlayer tunneling with a strength w.In H τ , r and k are, respectively, position and momentum operators, m * is the effective mass, the momentum shifts κ ± = [4π/(3a M )] (− √ 3/2, ∓1/2) are located at corners of the moiré Brillouin zone, and g j = 4π/( √ 3a M ) {cos[(j − 1)π/3], sin[(j − 1)π/3]} for j = 1, ..., 6 are the moiré reciprocal lattice vectors.Here a M ≈ a 0 /θ is the moiré period, and a 0 is the monolayer lattice constant.The moiré Hamiltonian H τ is constructed based on continuum approximation.Effects of in-plane lattice relaxation, which can be important at small twist angles (θ ≲ 1 • ), are not taken into account in the Hamiltonian H τ .The continuum approximation can also become less accurate at large twist angles (θ ≳ 5 • ) where the superlattices become small.Nevertheless, the Hamiltonian H τ qualitatively captures the evolution of electronic states as a function of θ in a convenient way.
The moiré Hamiltonian H τ with a fixed τ index is invariant under the C 3z and the combined C 2y T symmetry operations, where T is the time-reversal symmetry.The two valleys are related by the C 2y and the T symmetries, which map H + to H − .Moiré bands of H τ can be characterized by Chern numbers C + of the first band in +K valley is always quantized to 1.However, (C 1 , where θ * 1 ≈ 1.74 • .The topological phase transition at θ * 1 is accompanied by the gap closing between the second and the third bands at the γ point (i.e., the center of the moiré Brillouin zone).The first three moiré bands with a total Chern number of 0 provide a low-energy approximation of the full moiré band structure and are the focus of this work.

B. Wannier states
We construct Wannier states for the first three moiré bands so that the bands can be represented in terms of local orbitals.The center of the Wannier states can be determined by the C 3z eigenvalues of Bloch states at the high-symmetry momenta (i.e., γ and κ ± ).Symmetry analysis of the moiré Hamiltonian and Bloch states is performed in Appendix A following the approach developed in Ref. 44.The C 3z angular momenta of the Bloch states at γ and κ ± in the first three bands are labeled in Figs.1(b) and 1(c).The symmetry eigenvalues indicate that Wannier states for the first three bands should be centered at the A, B, and O positions, respectively.
We complement the symmetry analysis with a physical picture of the Wannier states.We first consider the case with θ < θ * 1 .In this first case, the first two bands have opposite Chern numbers (±1), and the third band has a Chern number of 0. Previous studies [2,11] have shown that the first two bands in one valley effectively realize a generalized Haldane model on a honeycomb lattice with two Wannier states located at the A and B positions, respectively.For the third band, the C 3z eigenvalues take the same value at γ and κ ± ( Fig. 1 (b)), which implies that the corresponding Wannier states are centered at the O positions.These considerations are consistent with the results obtained from the symmetry analysis.We then consider the second case with θ > θ * 1 , which is separated from the first case by a band inversion between the second and the third bands.This band inversion does not change the centers of Wannier states for the three bands as a whole.
We construct the Wannier states for θ < θ * 1 and θ > θ * 1 in a unified manner.The Wannier states in the moiré unit cell associated with the lattice site R = 0 can be formally expressed as, where τ is the valley index, n labels the three Wannier orbitals (n =1,2,3), and N is the number of k points included in the summation.Wannier states at a generic site R are obtained through lattice translation, In Eq. ( 2), ϕ nτ k (r) is a Bloch-like state associated with the Wannier states and related to Bloch states of H τ by a unitary transformation, where ψ n ′ τ k (r) is the Bloch state of the n ′ th moiré band at momentum k and U kτ is a 3×3 unitary matrix.We obtain U kτ by requiring that ϕ 1τ k and ϕ 3τ k are, respectively, maximally polarized to the top and bottom layers.This maximum value problem is equivalent to find the eigenstates of the layer polarization operator σ z (i.e., the z Pauli matrix in the layer pseudospin space) projected to the subspace spanned by ψ nτ k with n = 1, 2, 3, where T is a two-component spinor in the layer-pseudospin space.U kτ is then determined by, where We further fix the phase of ϕ nτ k in the following way.The phase of ϕ 1τ k is chosen by requiring that its top layer component ϕ 1τ kt is real and positive at r A = a M (1/ √ 3, 0), which is the A site associated with R = 0. Similarly, we take the bottom layer component of ϕ 3τ k to be real and positive at the B site r B = a M (2/ √ 3, 0).For ϕ 2τ k , we first fix the gauge at k = γ according to ϕ 2τ γb (r O ) > 0, where r O = a M (0, 0); then, at a generic k, we require that For convenience, hereafter we denote the Bloch-like states ϕ nτ k with n = 1, 2, and 3, respectively, as ϕ Aτ k , ϕ Oτ k , and ϕ Bτ k , and similarly for the Wannier states W nτ .
The Wannier states constructed for θ = 2.5 • using the above methods are shown in Fig. 2, where the amplitude of each layer component of W nτ is plotted.The Wannier states W Aτ and W Bτ are, respectively, polarized mainly to the top and bottom layer, and transformed to each other through the C 2y T operation (plus a lattice translation).The Wannier state W Oτ has significant weights on both layers and is an eigenstate of the C 2y T symmetry.
A tight-binding Hamiltonian based on the Wannier states can be formally constructed in real space and then Fourier transformed to obtain the momentum-space Bloch Hamiltonian [44], which has the following secondquantized form in the basis of Here c † knτ (c knτ ) is the electron creation (annihilation) operator of the Bloch-like state ϕ nτ k , and E nτ k is the energy of the Bloch state ψ nτ k with respect to H τ .The Hamiltonian Ĥ0 presents a noninteracting threeorbital model.It faithfully captures the symmetry and topology of the first three moiré bands.By diagonalizing h kτ , we reproduce the band structure of the first three bands and also obtain their weights projected onto the Wannier orbitals, as shown in Figs. 1 (b) and (c).For θ < θ * 1 , the first two moiré bands are mainly derived from the hybridization of A and B orbitals, while the third band is mainly composed of O orbital.For θ > θ * 1 , the second and the third bands exchange their orbital characters at the γ point, which results in the band inversion.

C. Interacting model
We now construct the interacting three-orbital model.For the charge-neutral twisted homobilayer, all valence band states are below the Fermi energy.When the system is doped with holes, it is more convenient to use the hole basis.We define the hole operator as b knτ = c † knτ .In the hole basis, the noninteracting Hamiltonian Ĥ0 is equivalent to where a constant term is dropped.The Coulomb interaction projected onto the three Wannier orbitals in the hole representation is expressed as where the summation is over the momentum k j (summed over the moiré Brillouin zone), the orbital index n j , and the valley index τ .The Coulomb matrix element is given by, where A is the system area.The structure factor where ℓ is the layer index, and φnτ k (r) = [ϕ nτ k (r)] * due to the particle-hole transformation.We use the dual-gate screened Coulomb interaction with the momentumdependent potential V q = 2πe 2 tanh (|q|d)/(ϵ|q|), where d is the gate-to-sample distance and ϵ is the dielectric constant.Here we neglect the effect of the vertical separation d 0 between the two MoTe 2 layers on the Coulomb potential since d 0 is much smaller than the typical interparticle distance.
By combining the noninteracting Hamiltonian Ĥ0 in Eq. ( 7) with the Coulomb interaction term Ĥint in Eq. ( 8), we obtain the full Hamiltonian of the interacting three-orbital model, Here we formulate the Hamiltonian Ĥ in momentum space.An alternative approach would be to construct the Hamiltonian in real space by projecting the Coulomb interaction onto the Wannier states.This latter approach works well to capture onsite Hubbard interactions and offsite density-density interactions [10], but is not suited to keep track of all relevant non-local interactions such as intersite Hund's coupling, pair hopping, and interactionassisted hopping terms [48][49][50].As shown in Fig. 2, neighboring Wannier states can have significant spatial overlap, and can lead to enhanced non-local interactions that play an important role in determining the ground state properties [48][49][50].Therefore, we do not use the real-space formalism.The momentum-space Hamiltonian Ĥ implicitly includes all microscopic interactions, and has been widely employed in moiré systems such as twisted bilayer graphene [51] to study correlated states [52,53].

III. PHASE DIAGRAM
We perform self-consistent mean-field studies of the Hamiltonian H using the Hartree-Fock (HF) approximation for integer filling factors ν = 1 and 2 at zero temperature.Here ν = 1 N k,n,τ b † knτ b knτ is the number of holes per moiré unit cell.At a given filling factor, we study multiple mean-field ansatzes that can break various symmetries of the system.The self-consistent calculation starting from different ansatz states often generates different mean-field solutions, whose energetic competitions determine the mean-field ground state.
We calculate the phase diagram as a function of the twist angle θ, the gate-to-sample distance d, and the dielectric constant ϵ.The parameter d can be experimentally varied by controlling the thickness of the encapsulating hBN layer.The screening effect by metallic gates is enhanced by reducing d.The dielectric constant ϵ accounts for the environmental screening from hBN as well as internal screening from remote moiré bands.Here we take ϵ as a phenomenological parameter and vary it theoretically to reveal the most generic phase diagram.nm and ϵ = 15, the ground state is in (1) the LP-FM z phase for θ < 2 • , (2) the quantum anomalous Hall insulator (QAHI) phase for 2 • < θ < 2.9 • , (3) the topologically trivial FM z phase for 2.9 • < θ < 3.9 • , and (4) the O − 120 • AF ± phase for 3.9 • < θ < 5.0 • .Here the LP-FM z phase is a multiferroic state with coexisting ferroelectricity and ferromagnetism.The ferroelectricity arises from spontaneous layer polarization (LP) induced by imbalanced hole densities on A and B Wannier orbitals [red line in Fig. 4(d)], which is driven by the Coulomb repulsion between A and B sites.The magnetism is a result of valley polarization, which leads to out-of-plane ferromagnetism (FM z ) since valleys are locked to out-of-plane spins in the system.We note that the LP-FM z phase has been previously obtained in an exact-diagonalization study of tMoTe 2 for θ limited to be below 1.5 • [16], and multiferroic states have also been studied in moiréless graphene systems [54,55].The next two phases in our phase diagram, QAHI and FM z , are both valley polarized but layer unpolarized; they are distinguished by the total Chern number C of the ground state, which is 1 in QAHI and 0 in FM z .The QAHI and FM z phases are separated by a weakly first-order phase transition from our numerical calculations, which is manifested by a slight jump of the occupation number of O orbitals across the transition [Fig.4(d)].Since this transition is first order, the gap does not need to completely close at the transition, as shown in Fig. 4(d).The topological transition is also revealed by the HF band struc-ture in the two phases, as shown in Figs.4(b) and 4(c).The band (in the original basis defined by c † knτ and c knτ operators) above the Fermi level is derived from A and B orbitals at γ and κ ± momenta in the QAHI phase, but is composed of O orbital at γ point in the FM z phase.This band inversion results in the change of the Chern number.Lastly, the O − 120 • AF ± phase has the in-plane 120 • antiferromagnetic (AF) Néel order on the O sites, and the superscript ± indicates the vector chirality of the antiferromagnetism [9], as illustrated in Fig. 3c.This AF state with a given vector chirality spontaneously breaks the C 2y symmetry, but states with opposite chiralities are related by the C 2y T symmetry and therefore, energetically degenerate.Because of the C 2y symmetry breaking, the O − 120 • AF ± states also have spontaneous layer polarization controlled by the vector chirality of the antiferromagnetism.This magnetoelectric coupling allows electric tuning of the vector chirality.The layer polarization is numerically confirmed by the imbalanced densities on A and B sites in the O − 120 • AF ± phase, as demonstrated by the red line in Fig. 4(d).
The phase diagram is determined by comparing the energies of multiple competing states, as shown in the lower panel of Figs.3(a) and 3(b).In addition to the states described above, we also find LP-FM x , LP-120 • AF ± , and O-FM x states as mean-field solutions in certain parameter regimes, although these latter states are less energetically favorable and do not appear in the ground-state phase diagram.Here LP-FM x refers to a state with spontaneous layer polarization and in-plane ferromagnetism (FM x ) on either A or B sites.LP-120 • AF ± is also layer polarized, but has in-plane 120 • AF order on either A The phase transition between competing phases is generically first order in the ν = 1 phase diagram.Across the transitions, the energy of the ground state is continuous, but other quantities, such as the sublattice-resolved hole density, layer polarization, and gap, can have discontinuous jumps.The evolution of these quantities as a function of θ for the ground state is shown in Fig. 4(d).The θ dependence of the phase diagram reflects the competition between the kinetic energy and the Coulomb interaction.As the moiré period a M decreases by increasing θ, the ratio of the characteristic Coulomb repulsion e 2 /ϵa M to the single-particle bandwidth typically decreases, and different phases can appear.
We turn to the ϵ = 10 phase diagram in Fig. 3(a), which includes only the LP-FM z phase for θ ≲ 2.4 • and the FM z phase for θ ≳ 2.4 • .The QAHI state is a solution to the HF equations for θ ≲ 2 • , but it is energetically unfavorable compared to the LP-FM z phase in this small angle regime.The stronger Coulomb interaction for ϵ = 10 tends to localize carriers by repulsion.The holes are localized at A or B sites in the LP-FM z phase, but around O sites in the FM z phase.The existence of the QAHI phase requires an intermediate strength of the Coulomb repulsion.
The evolution of the phase diagram as a function of ϵ is shown in Fig. 5.The region of the FM z phase in the parameter space of (θ, ϵ) shrinks as ϵ increases and finally vanishes for ϵ > 18, where there is a first-order phase transition between the QAHI phase and the O − 120 • AF ± phase.The phase boundary between the QAHI phase and the LP-FM z (O − 120 • AF ± ) phase shifts to lower θ as ϵ increases.The reason is that the interaction effects become weaker with increasing ϵ but stronger with decreasing θ.The O − 120 • AF ± phase can become gapless at very large ϵ and finally give way to the metallic state without any symmetry breaking.
We now discuss the d dependence of the phase diagrams.The screening from the gates mainly suppresses the Coulomb repulsion for the interaction range larger than d.When d ≫ a M , the system should have a weak dependence on d since the gates only screen the longrange tail of the Coulomb repulsion in this parameter regime.Numerically, we find that the phase boundaries indeed have a relatively weak d dependence for d down to 10 nm, as shown in the phase diagrams of Figs.3(a   We present additional discussion on the ν = 1 phase diagram in Appendices B and C. In Appendix B, we show that (1) the QAHI state is captured by the Kane-Mele-Hubbard model, where the Hubbard interaction combined with the spin-orbit coupling drives the out-of-plane spin polarization at ν = 1; (2) the transition from the LP-FM z to the QAHI state is phenomenologically described by a spinless Haldane model with nearest-neighbor repulsion; (3) the O − 120 • AF ± state can be understood as a result of Fermi surface instability.These analyses deepen the physical understanding of the phase diagram.In Appendix C, we compare results obtained, respectively, from one-band, two-band, and three-orbital models, which show that the three-orbital model is essential to capture all the competing states.(3) the spin Hall insulator (QSHI) phase for θ 4.5 • .The phase diagram is obtained by comparing the energy of multiple competing states, including FM z , FM x , AF x , z , and symmetric states.All these states equal hole densities on A and B sublattices.The state does not break any symmetry of the system and always realizes the QSHI state for the θ range (1.5 • , 5 • ) consideration.In the z (FM x ) state, A and B sublattices have ferromagnetic moments along z (x) direction.The AF z (AF x ) has collinear antiferromagnetic Néel order with opposite magnetic moments on A and B sublattices along z (x) direction.The HF band structures for the FM x , AF x , and phases are, respectively, shown in Figs.7(a In the small θ of θ < 1.7 • , the the FM z state is than that of AF x state only by a very small margin, as shown by the bottom panel of Fig. 6(b).The transition between FM z and AF x phases at θ ≈ 1.7 • is order.By contrast, the tranbetween the AF x and phases θ 4.5 is second order.This is manifested by the continuous supof the Néel order as θ approaches 4.5 • from AF x phase [blue line in Fig. 7(d)].This phase transition between the AF x and QSHI phases is reminiscent of the same transition in the Kane-Mele-Hubbard model at a critical onsite Hubbard interaction [56].While our realistic model differs from the Kane-Mele-Hubbard model due to complications such as the additional O orbital and the long-range Coulomb interaction, the latter model provides important intuitions in understanding the phase diagram.Mean-field theory underestimates the critical Hubbard repulsion in the Kane-Mele-Hubbard model compared to that obtained from quantum Monte Carlo simulation [57].It is likely that our mean-field phase diagram also overestimates the tendency towards the formation of the Néel order, and the QSHI phase could appear in a larger parameter space.
In the QSHI phase, the occupied bands in opposite valleys have opposite Chern numbers of ±1.Since the two valleys also act as spin up and down that are related by the T symmetry, the QSHI state is characterized by nontrivial Z 2 topological invariant.In the single-particle band structure of the noninteracting moiré Hamiltonian, the first and the second band overlap in energy for θ > 3.1 • , and the system would then be a metallic state at ν = 2 without taking into account of interaction effects.However, the QSHI phase of the interacting model has a full charge gap as indicated by the HF band structure shown in Fig. 7(c).Thus, the Coulomb repulsion is crucial in opening up the charge gap in the QSHI phase.
When ϵ is reduced to 10, the corresponding phase diagram in Fig. 6(a) only has the FM z and AF x phases for θ ∈ (1.5 • , 5 • ).The stronger Coulomb interaction for ϵ = 10 enhances the antiferromagnetic order and therefore, the QSHI state no longer appears in the phase diagram of Fig. 6(a).The existence of the QSHI state requires an intermediate interaction strength, which should be strong enough to open a full charge gap but weak enough to avoid the instability toward magnetism.This is clearly shown by the phase diagram as a function of θ and ϵ plotted in Fig. 8.

IV. DISCUSSION AND CONCLUSION
In summary, we present a systematic study of the quantum phase diagrams of tMoTe 2 that result from the interplay between band topology and many-body interaction.At the single-particle level, we construct a threeorbital model that can be viewed as a generalization of the two-orbital Kane-Mele model.The three-orbital model faithfully captures the band dispersion, symmetry, and topology of the low-energy moiré bands for a large range of θ.For the interacting physics, we obtain the quantum phase diagrams based on comprehensive manybody calculations by comparing the energy of multiple competing states.At ν = 1, we find a cascade of phase transitions as θ increases.In the small θ regime, the system is in a layer-polarized phase with spontaneous imbalanced hole densities on A and B sites.The QAHI phase appears in an intermediate range of θ and is superseded by topologically trivial phases (i.e., FM z and O − 120 • AF ± phases) in the large θ regime.At ν = 2, the QSHI phase can appear in the large θ regime, but undergo an instability towards an antiferromagnetic state below a critical θ.
The obtained phase diagrams serve as a guide for further experimental and theoretical studies, but the results should be taken to be qualitative instead of quantitative.There are several theoretical difficulties in establishing truly quantitative phase diagrams.First, parameters in the moiré Hamiltonian are subjected to numerical uncertainties.We take the parameter values from Ref. [2], where parameters were obtained by fitting to the density-functional-theory (DFT) band structures of the untwisted homobilayer at different high-symmetry stackings.The DFT calculations are likely not accurate enough to pin down the exact values of small quantities such as V and w that are on the order of 10 meV.Second, mean-field theory can capture the quantum phase diagram qualitatively but may overestimate the tendency towards symmetry-breaking states.In spite of these quantitative issues that are generic challenges for moiré systems, we expect that our study captures the qualitative features of the phase diagrams.
Recent experiments on tMoTe 2 reported optical signatures of a ferromagnetic correlated insulator in a θ = 3.9 • sample [21] and a quantum anomalous Hall state with Chern number of 1 in another θ = 3.7 • sample [22], where the hole filling factor ν is 1 in both cases.The topological character of the ν = 1 ferromagnetic insulator in the θ = 3.9 • sample was not reported [21].An independent experiment reported the quantum anomalous Hall state in tMoTe 2 with θ = 3.4 • at ν = 1 [23].Our ν = 1 phase diagram featuring both the topological ferromagnetic phase (i.e., the QAHI phase ) and the topologically trivial ferromagnetic phase (i.e., the FM z phase) is qualitatively consistent with these experimental findings.In the phase diagram of Fig. 3(b), the QAHI phase appears for θ in the range of about (2 • , 2.9 • ), which is below the experimental twist angle 3.4 • .This discrepancy can be remedied by slightly adjusting the parameters used in the moiré Hamiltonian.In Fig. 9, we present the ν = 1 phase diagram as a function of θ and V , where V is the moiré potential amplitude as defined in Eq. ( 1).In the calculation of Fig. 9, we keep other parameters (except θ and V ) of the moiré Hamiltonian in Eq. ( 1) fixed, and take ϵ = 15 and d = 20 nm.As shown in Fig. 9, the θ range for the QAHI phase is shifted to (2.35 • , 4.05 • ) if V is increased to 9 meV, which leads to better consistency with the experimental finding [22].Increasing the intralayer moiré potential amplitude V has two effects.First, a deeper moiré potential leads to stronger spatial confinement of A and B orbitals, and therefore, weaker hopping parameters between A and B sites.The doped holes then have a stronger tendency to be localized by Coulomb repulsion, which shifts the phase boundary between the LP-FM z and QAHI phase to larger twist angles.Second, the deeper intralayer moiré potential also increases the single-particle energy difference between A/B orbital and O orbital, and effectively pushes the O orbital away from the Fermi level.Therefore, the topologically trivial FM z phase gradually shrinks as V increases and finally disappears at V ≈ 8.7 meV.For V > 8.7 meV, there is a first-order phase transition between the QAHI phase and the O − • AF ± phase, where the corresponding boundary also shifts towards larger twist angles as V increases.We note that moiré potential could in principle be tuned by pressure [58,59], but the exact dependence of model parameters and the phase diagrams on the pressure requires future investigations.
We comment on the implication of our results for tWSe 2 .Early transport studies of tWSe 2 reported topologically trivial correlated insulator at ν = 1 in the twist angle range from 4 • to 5.5 • [18,19].Recent local electronic compressibility measurements of tWSe 2 at ν = 1 revealed QAHI state in a narrow twist angle range between 1.2 • and 1.25 • , but topologically trivial correlated insulator for θ either below or above this range [20].Qualitatively, our ν = 1 phase diagram with the QAHI phase appearing in an intermediate range of θ is consistent with these experimental observations.The O − 120 • AF ± phase and the LP phase can capture the topologically trivial correlated insulators observed at large and small twist angles, respectively.Quantitatively, one peculiar experimental finding is that the QAHI state in tWSe 2 only appears at a very narrow twist angle range [20], which calls for further investigation.
Signatures of fractional quantum anomalous Hall (FQAH) states at fractional filling factors were also reported in tMoTe 2 [22,23].Exact diagonalization studies of the FQAH states have recently been performed, where the calculations were done in the truncated Hilbert space projected to the first moiré valence band [60,61].Our study highlights the importance of including multiple bands in the many-body problem.For example, the fully valley polarized state in the truncated Hilbert space projected to the first moiré valence bands is always a quantum anomalous Hall state at ν = 1, but this is not the case in the three-orbital model (see more discussion in Appendix C).It is desirable to examine the energy competition between the FQAH state and other candidate states (e.g., charge density wave state) at fractional filling factors in the multi-band model using techniques beyond the mean-field theory.
FQAH states are analogous to fractional quantum Hall states formed in Landau levels.However, there is a crucial difference between Landau levels and topological moiré bands in tMoTe 2 , since the latter respects timereversal symmetry.The FQAH states realized in tMoTe 2 spontaneously break the time-reversal symmetry.Timereversal symmetric fractional quantum spin Hall insulators [62,63], as another distinct type of fractionalized states, may also emerge in tMoTe 2 and tWSe 2 .
The obtained quantum phases can have interesting responses to external fields.For example, an out-of-plane electric displacement field can drive the QAHI state at ν = 1 into a topologically trivial layer polarized state.An out-of-plane magnetic field turns the ν = 1 O − 120 • AF ± state into a canted antiferromagnetic state, and finally to a ferromagnetic state for magnetic field above a critical value, where the ferromagnetic state can be in the QAHI phase, as shown in Fig.The exploration of the rich quantum states in tMoTe 2 is just at its beginning stage.Note added.Fractionally quantized anomalous Hall effect has recently been observed in tMoTe 2 through transport measurement [64,65].
Appendix A: C3z symmetry eigenvalues The single-particle wave function can be viewed as a direct product of the spatial envelope wave function ψ nτ k (r), the atomic state |d τ ⟩, and the spin state |s τ ⟩, where |d τ ⟩ = |d x 2 −y 2 ⟩ + iτ |d xy ⟩ stands for the dominant atomic orbital, and |s τ ⟩ represents spin up and down, respectively, for τ = + and −.The C 3z symmetry acts separately on the three parts.For the internal degrees of freedom, the C 3z symmetry eigenvalue is given by, where l τ = −τ /2 is the angular momentum contributed by the spin and atomic orbital.
We now study the C 3z symmetry property of the spatial wave function using the approach developed in Ref. 44.We first apply the following unitary transformation to H τ to make the C 3z symmetry transparent, where ∆T,τ (r) = w(e iτ q1•r + e iτ q2•r + e iτ q3•r ).Here q 1 = κ − − κ + , q 2 = R3 q 1 , q 3 = R3 q 2 , and R3 is the real-space anticlockwise rotation by 2π/3.The new Hamiltonian Hτ (r) is clearly invariant under the threefold rotation of r since Hτ ( R3 r) = Hτ (r).
The C 3z symmetry of the original Hamiltonian H τ can be represented as Ĉ3z,τ = Λ −1 τ (r) R3 Λ τ (r), which leads to the symmetry identity Ĉ3z,τ H τ (r) Ĉ−1 3z,τ = H τ (r).The C 3z symmetry acts on the wave function in the following way, where ψnτ k (r) = Λ τ (r)ψ nτ k (r).At the three high-symmetry momenta γ and κ ± , the Bloch wave function is an eigenstate of the C 3z symmetry, where L nτ k is the C 3z angular momentum of the spatial wave function and defined modulo 3. Using Eqs.(A3) and (A4), we find that Therefore, L nτ k can be determined by numerically comparing ψnτ k ( R3 r) and ψnτ k (r).The values of L nτ k for the first three bands at the three high-symmetry momenta are labeled in Figs.1(b) and 1(c).
In a system with C 3z symmetry [66], the Chern number of a band is determined by the sum of the C 3z angular momentum up to modulo 3, τ ] mod 3 = 0.The calculated Chern numbers and C 3z angular momentum are consistent with this condition.
B2) where the lattice translational symmetry is assumed to be preserved, N is the number of unit cells, and ℓ is the sublattice index.
By diagonalizing the mean-field Hamiltonian, we obtain the total energy at ν = 1 as follows, where and M ℓ is parameterized as M ℓ (sin θ ℓ cos ϕ ℓ , sin θ ℓ sin ϕ ℓ , cos θ ℓ ).In deriving the total energy, we assume that n A = n B and M A = M B = M .It is obvious that E tot is minimized when ϕ A = ϕ B .By further taking ∂E tot /∂θ ℓ = 0, we find two types of solutions, (1) Therefore, the Kane-Mele-Hubbard model supports the out-of-plane spin polarized state at ν = 1, which exactly gives rise to the valley polarized QAHI state in tMoTe 2 .From the expression of E tot , we see that the spin-dependent term F z σ z s z in the Kane-Mele model leads to Ising anisotropy and stabilizes the QAHI state.

Haldane-V model
The spontaneous sublattice (layer) polarization in the LP-FM z phase is driven primarily by the repulsion between A and B orbitals.To capture this physics, we focus on states with out-of-plane spin polarization.In this spin sector, we can study the spinless Haldane model with nearest-neighbor repulsion V , which is expressed as, where ĤH is the Haldane model on the honeycomb lattice and ĤV describes the nearest-neighbor repulsion.We note that ĤKM can be understood as two copies of ĤH that are time-reversal partners.
In the Hartree-Fock theory, Ĥ can be approximated as, where b k = {b kA , b kB } T and f ∥ = F ∥ /t 0 .We define a vector ρ as ⟨b + σb⟩ /2, where b = {b A , b B } T for a pair of nearest neighbors.ρ z is the z component of ρ and quantifies the sublattice polarization.ρ ∥ is defined as ρ x − iρ y and quantifies the inter-sublattice coherence.The nearest-neighbor hopping parameter is renormalized from the bare value t 0 to be t 0 − V ρ ∥ by the interaction.
The total mean-field energy at ν = 1 is, We obtain the following two selfconsistent equations by, respectfully, taking ∂Etot ∂ρz = 0 and ∂Etot ∂ρ ∥ = 0, We note that the self-consistent equations always have a symmetric solution with no sublattice polarization (ρ z = 0), but it can also have a symmetry-breaking solution with spontaneous sublattice polarization (ρ z ̸ = 0).By minimizing the total energy, we obtain the phase diagram [Fig.10] of the Haldane-V model as a function of V /t 0 and t ′ 1 /t 0 , where t ′ 1 is the imaginary part of one of the second nearest neighbor hopping parameters.For a fixed t ′ 1 /t 0 , ρ z is 0 (nonzero) for V below (above) a critical value.Therefore, spontaneous sublattice polarization can be driven by V above a critical value.This qualitatively explains why the LP-FM z becomes the ground state at small twist angles, as V /t 0 increases with decreasing twist angle.The state with a spontaneous finite ρ z can also be described as a charge density wave state with imbalanced charge density at A and B sites.In tMoTe 2 , A and B orbitals have opposite layer polarization, and therefore, spontaneous sublattice polarization also leads to spontaneous layer polarization.

Fermi surface instability
At large twist angles, the Coulomb interaction can become less dominant compared to the single-particle kinetic energy, and the fermiology associated with the Fermi surface can become important.In Fig. 11, we plot the noninteracting Fermi surface at half-filling (i.e., ν = 1) of the topmost moiré valence bands for θ = 4.5 • .There is an approximate Fermi surface nesting, and the nesting vector Q is approximately equal to κ − − κ + , which can drive Fermi surface instability with a √ 3 × √ 3 real-space texture.This instability exactly explains the formation of the O −120 • AF ± state at large twist angles.The QAHI state is the same ground state for one-band, two-band model, and three-orbital models for a certain range of twist angles.However, there is a crucial difference in the size of the charge gap.As shown in Fig. 13, the one-band and two-band models generally overestimate the charge gap compared to the three-orbital model.In the three-orbital model, the charge gap has a chargetransfer nature, since it is determined by electron and hole located at different real-space orbitals but within the same valley.The one-band model completely fails to capture this charge-transfer gap, and therefore, overestimates the charge gap.
We also perform the Hartree-Fock calculation in the plane-wave basis without projecting onto the first three bands, which leads to quantum phase diagrams that quantitatively agree with that from the three-orbital model, as we report in a follow-up work [67].Therefore, the three-orbital model can be viewed as a minimal model that provides a unified description of competing phases over a large range of twist angles.

FIG. 1 .
FIG. 1.(a) Moiré superlattices formed in tMoTe2.The dots with green, red, and blue colors indicate, respectively, O, A, and B sites, which also correspond to the centers of the three Wannier orbitals.The black dashed lines mark a moiré unit cell.(b), (c) Moiré band structure at θ = 1.5 • and 2.5 • in +K valley.The integer numbers at γ and κ± label the C3z angular momentum of the first three bands (see Appendix A), while the color represents the relative weight of the three Wannier orbitals.The band structures are plotted along highsymmetry paths in momentum space with κ ′ + = −κ−.
n is the band index (counting in descending order of energy).Because of the T symmetry, C on the single-particle topological phase diagram of H τ as a function of the model parameters (V, ψ, w) can be found in Ref. 9. Unless otherwise stated, we take the parameter values estimated for tMoTe 2 from Ref. 2, a 0 = 3.472 Å, m * = 0.62m e , V = 8 meV, ψ = −89.6 • , and w = −8.5 meV, where m e is the electron bare mass.The calculated moiré band structures are plotted in Figs.1(b) and 1(c) for two representative twist angles θ = 1.5 • and 2.5 • , respectively.For twist angles θ ≤ 5 • studied in this work, C

FIG. 3 .
FIG. 3. (a), (b) Top panel: The quantum phase diagram at ν = 1 as a function of θ and d.The color map plots the charge gap.White solid lines mark the first-order phase transitions.Bottom panel: Energy per moiré unit cell of competing states relative to the layer unpolarized but valley polarized (LUP-VP) states at d = 20 nm.ϵ is 10 in (a) and 15 in (b).(c) Schematic illustration of competing states.In the LP-FMx, LP-FMz, LP-120 • AF + and LP-120 • AF − states, the A and B sublattices have unequal hole densities.In LP-FMx (LP-FMz) state, A (or B) sublattices have ferromagnetism along x (z) direction.In LP-120 • AF + (LP-120 • AF − ) state, the Néel order has an anticlockwise (clockwise) arrangement on A (or B) sublattices along the y direction.Similar notations are used for O-FMx and O-120 • AF ± states, of which magnetic moments are developed on the O sublattice.In QAHI and FMz states, the A and B sublattices have equal density and the magnetization is along the z direction.The total Chern number is 1 in QAHI and 0 in FMz states, both of which can be classified as the LUP-VP states.

FIG. 4 .
FIG. 4. (a)-(c) The HF band structure at ν = 1 in (a) the LP-FMz phase at θ = 1.5 • , (b) the QAHI phase at θ = 2.5 • , and (c) the FMz phase at θ = 3.5 • .The band structures are presented in the basis defined by c † knτ and c knτ operators.The solid and dotted lines, respectively, plot bands from +K and −K valleys.The color represents the weight of the three Wannier orbitals.The middle of the interaction-induced gap, marked by the black dashed line, is set to 0 in each plot.The integer numbers at γ and κ± label the C3z angular momentum of the band above the Fermi energy.(d) The sublattice polarization PAB (red line), the hole number nO on each O site (blue dashed line), and the charge gap (black line) in ν = 1 ground states as a function of θ.PAB characterizes the layer polarization and is defined as |nA − nB|/(nA + nB), where nA (nB) is the hole number on each A (B) site.The vertical dashed lines mark the transitions between different phases.d is 20 nm and ϵ is 15 for (a)-(d).

FIG. 5 .
FIG. 5.The quantum phase diagram at ν = 1 and d = 20 nm as a function of θ and ϵ.The color map plots the charge gap.The symbol "M" stands for the metallic state without any symmetry breaking.The white dashed line separates gapped and gapless O-120 • AF ± states.
) and 3(b).For small d, the phase boundary between the LP-FM z and QAHI phases in the ϵ = 15 phase diagram of Fig.3(b) shifts towards the LP-FM z side as d decreases from 10 nm to 5 nm; the Coulomb repulsion between the nearest A and B sites is reduced under such small values of d, and the LP-FM z phase becomes less energetically favorable.A similar effect is found for the phase boundary between the FM z and QAHI phases, as shown in

FIG. 6 .
FIG. 6. (a), (b) Top panel: The quantum phase diagram at ν = 2 as a function of θ and d.The color map plots the charge gap.White solid (dashed) lines mark the first (second)-order phase transitions.Bottom panel: Energy per moiré unit cell of competing states relative to the symmetric (sym) state at d = 20 nm.ϵ is 10 in (a) and 15 in (b).(c) Schematic illustration of competing states at ν = 2.

Fig. 3 (
Fig. 3(b).The numerical results indicate that the parameter regimes of the QAHI phase can be slightly enlarged by reducing d down to below a M .We present additional discussion on the ν = 1 phase diagram in Appendices B and C. In Appendix B, we show that (1) the QAHI state is captured by the Kane-Mele-Hubbard model, where the Hubbard interaction combined with the spin-orbit coupling drives the out-of-plane spin polarization at ν = 1; (2) the transition from the LP-FM z to the QAHI state is phenomenologically described by a spinless Haldane model with nearest-neighbor repulsion; (3) the O − 120 • AF ± state can be understood as a result of Fermi surface instability.These analyses deepen the physical understanding of the phase diagram.In Appendix C, we compare results obtained, respectively, from one-band, two-band, and three-orbital models, which show that the three-orbital model is essential to capture all the competing states.

B. ν = 2
phase diagram The quantum phase diagram at ν = 2 as a function of θ and d is shown in Figs.6(a) and 6(b) for ϵ = 10 and 15, respectively.The phase diagram again has a weak d dependence.Therefore, we focus on the θ dependence of the phase diagram.The ground state along line cut with d = nm ϵ = is (1) the FM z phase for θ < 1.7 • , (2) the AF x for 1.7 • < θ < 4.5 • , ), 7(b), and 7(c), where the color encodes the weight on the O orbital.results show that the holes at ν = 2 are primarily doped to the A and B sites.

FIG. 7 .
FIG. 7. (a)-(c) The HF band structure at ν = 2 in (a) the FMz phase at θ = 1.5 • , (b) the AFx phase at θ = 2.5 • , and (c) the QSHI phase at θ = 4.7 • .The band structures are presented in the basis defined by c † knτ and c knτ operators.The color represents the weight of the O orbital.In (a), the solid and dotted lines, respectively, plot bands from +K and −K valleys.In (b) and (c), each band is doubly degenerate.(d) The charge gap (black line) and the in-plane magnetic moment M A x on A sublattice (blue line) in ν = 2 ground states as a function of θ. d is 20 nm and ϵ is 15 for (a)-(d).

FIG. 8 .
FIG. 8.The quantum phase diagram at ν = 2 and d = 20 nm as a function of θ and ϵ.The color map plots the charge gap.The symbol "M" stands for metallic state.
This work is supported by National Key Research and Development Program of China (Grants No. 2021YFA1401300 and No. 2022YFA1402401), National Natural Science Foundation of China (Grant No. 12274333), and start-up funding of Wuhan University.W.-X. Q. is also supported by Postdoctoral Innovation Research Position of Hubei Province (Grant No. 211000128).The numerical calculations in this paper have been done on the supercomputing system in the Supercomputing Center of Wuhan University. ℓ

E
FIG. 11.Noninteracting Fermi surface (white dashed line) at half-filling of the topmost moiré valence bands.The color map shows the energy dispersion in the first moiré Brillouin zone.The vector Q indicates an approximate nesting.

FIG. 12 .
FIG. 12. Energy per moiré unit cell of competing states relative to the LUP-VP states at ν = 1, ϵ = 15 and d = 20 nm.The results are obtained using (a) one-band and (b) two-band models, respectively.Other parameter values are the same as those used in Fig. 3(b).

FIG. 13
FIG. 13. (a)-(c)The HF band structure for the ν = 1 QAHI phase at θ = 2.5 • calculated with one-band, two-band, and three-orbital model, respectively.The solid and dotted lines, respectively, plot bands from +K and −K valleys.(d) The charge gap of the QAHI phase obtained from one-band (red), two-band (green) and three-orbital (blue) model for the twist angle range from 2 • to 2.9 • , where the QAHI phase is the same ground state in all three models.Parameter values are the same as those used in Fig.12.
Appendix C: Comparison of different modelsIn this Appendix, we present a comparison of results obtained, respectively, from one-band, two-band, and three-orbital models, as summarized in Fig.12.First, in the one-band model [Fig.12(a)],we project the interacting Hamiltonian onto the topmost moiré valence bands, where each valley contributes one band.As shown in Fig.12(a), the valley polarized state at ν = 1 is always the QAHI in the one-band model.This model fails to capture the LP-FM z state, the topologically trivial FM z state, and the O − 120 • AF ± state.(2)In the twoband model, we project the interacting Hamiltonian onto the first two moiré valence bands, where each valley contributes two bands.As shown in Fig.12(b), the twoband model captures the LP-FM z state for small twist angles and its transition to the QAHI state above a critical angle, but fails to describe the topologically trivial FM z and O − 120 • AF ± states.Quantitatively, the transition angle between the LP-FM z and QAHI states is underestimated in the two-band model compared to the three-orbital model shown in Fig.3(b).The additional O orbital in the three-orbital model is crucial in correctly finding all competing states such the topologically trivial FM z and the O − 120 • AF ± states.