Band topology, Hubbard model, Heisenberg model, and Dzyaloshinskii-Moriya interaction in twisted bilayer WSe$_2$

We present a theoretical study of single-particle and many-body properties of twisted bilayer WSe$_2$. For single-particle physics, we calculate the band topological phase diagram and electron local density of states (LDOS), which are found to be correlated. By comparing our theoretical LDOS with those measured by scanning tunneling microscopy, we conclude that the first moir\'e valence band is likely to be topologically trivial. For many-body physics, we construct a generalized Hubbard model on a triangular lattice based on the calculated single-particle moir\'e bands. We show that a layer potential difference, arising, for example, from an applied electric field, can drastically change the non-interacting moir\'e bands, tune the spin-orbit coupling in the Hubbard model, control the charge excitation gap of the Mott insulator at half filling, and generate an effective Dzyaloshinskii-Moriya interaction in the effective Heisenberg model for the Mott insulator. Our theoretical results agree with transport experiments on the same system in several key aspects, and establish twisted bilayer WSe$_2$ as a highly tunable system for studying and simulating strongly correlated phenomena in the Hubbard model.


I. INTRODUCTION
Twisted bilayers with a long-range moiré pattern provide highly tunable platforms to study fundamental physics for both single-particle and many-body phenomena.An important breakthrough was the experimental discovery of superconducting and correlated insulating states [1,2] in magic-angle twisted bilayer graphene (TBG) [3].While magic angle TBG is under active study and hosts a rich variety of phenomena [4][5][6][7], it poses challenges for both experiment and theory.In experiment, superconducting and correlated insulating states in TBG are fragile and appear only within a narrow range of twist angle around the magic angle (∼ 1.1 • ), requiring great experimental efforts to fine tune the twist angle.In theory, the low-energy moiré bands in TBG defies the construction of fully symmetric Wannier states because of intrinsic obstructions [8], which complicates theoretical analysis.
It was theoretically proposed that twisted bilayer transition metal dichalcogenides (TMDs) represent a simpler system compared to TBG and can provide a platform to simulate model Hamiltonians such as Hubbard model and Kane-Mele model [9,10].Here TMDs refer to group-VI semiconducting transition metal dichalcogenides such as WSe 2 [11].The simplicity of TMDs compared to graphene originates from the fact that the former is a semiconductor with a large band gap as well as a large spin-orbit coupling, while the latter is a semimetal with Dirac cones and spin SU(2) symmetry.Because of the reduced symmetries in TMDs, the low-energy degrees of freedom in twisted bilayer TMDs are fewer than TBG, which leads to theoretical simplification, allowing effective realizations of simple yet important model Hamil- * wufcheng@umd.edutonians [9,10].Another noticeable difference between twisted bilayer TMD and TBG is that the nearly flat moiré bands appear in a large range of twist angles in the former system, but only occur within a small window (±0.1 • ) around the magic angle in the latter system.This difference could lead to practical advantages, as there is no longer an acute need to carefully fine tune the twist angle in order to achieve the flat band situation.Singleparticle flat bands strongly enhance the relative interaction strength since the non-interacting kinetic energy is suppressed under the flat band condition, potentially leading to many interesting correlated quantum phases.
There are two types of twisted TMD bilayers, namely, heterobilayers and homobilayers.In heterobilayers, the two layers are respectively two different TMD materials, for example WSe 2 /MoSe 2 , which automatically lift the layer degeneracy.This moiré system can realize a generalized Hubbard model on a triangular lattice formed by effective moiré sites [9].Such a Hubbard model simulator based on TMD heterobilayers has recently been experimentally realized in Refs.12 and 13, which report evidence for Mott insulators and Wigner crystals.
In this paper we focus on twisted TMD homobilayers, where the two layers are formed from the same material.Because of stronger interlayer coupling, homobilayers can potentially be more interesting as well as more tunable compared to heterobilayers.Our work is motivated by two experimental studies on twisted bilayer WSe 2 , where one experiment is based on scanning tunneling microscope (STM) [14], and the other is on transport measurement [15].Both experimental papers [14,15] report signatures of narrow moiré bands in twisted bilayer WSe 2 , and the transport experiment [15] also identifies halffilled correlated insulators that can be sensitively tuned using an external displacement field.
The purpose of this work is mainly twofold.First, we determine the nature of the low-energy non-interacting moiré bands, including their topological character and their mapping to effective lattice models.We present systematic topological phase diagrams characterized by valley Chern numbers as a function of system parameters.We find that the topology of the first moiré valence band is closely connected with the pattern of electron density distribution in moiré superlattices.By comparing our theoretical local density of states with those measured by STM [14], we conclude that the first moiré valence band in twisted bilayer WSe 2 is likely to be topologically trivial, and can be described by a one-orbital tight-binding model on a triangular lattice.The tight-binding model combined with Coulomb repulsion leads to the realization of an effective Hubbard model for the corresponding interacting system.Second, we demonstrate the convenient tunability provided by an external out-of-plane displacement field in controlling both single-particle as well as many-body properties of TMD homobilayers.For singleparticle physics, we show that V z , a layer potential difference generated by the displacement field, drastically changes the moiré band structure, tunes van Hove singularities, and controls the effective spin-orbit coupling in the tight-binding model.For many-body physics, we predict that V z generates an effective Dzyaloshinskii-Moriya (DM) interaction in the effective Heisenberg model (the spin model for Mott insulator at half filling) associated with the Hubbard model, and acts as a tunable experimental knob that can turn on and off the corresponding correlated (Mott) insulators at half filling.Our theoretical results are consistent with a recent transport experiment in WSe 2 [15].
We highlight two specific important predictions of our theory.(1) While the first moiré valence band is likely to be topologically trivial, other moiré bands can still be topologically nontrivial.This should motivate transport study on the (topologically nontrivial) second and even third moiré valence bands by increasing the hole carrier density.(2) The DM interaction leads to vector spin chirality in the 120 • antiferromagnetic ground state of the Heisenberg model on a triangular lattice.This field-tunable DM interaction in the moiré system is an interesting phenomenon, which may find applications in spintronics.
The remainder of this paper is organized as follows.In Sec.II, we present a thorough study of moiré band structure in twisted bilayer WSe 2 with a focus on the topological character and the electron density distribution in real space.In Sec.III, we construct a tight-binding model for the first moiré valence band in the topologically trivial regime and in the presence of a finite V z .In Sec.IV, we construct a Hubbard model for the first moiré band by including Coulomb repulsion.We study the Hubbard model at half filling by mapping it to the corresponding Heisenberg model as well as directly by using a meanfield theory.The effects of V z on many-body physics are also calculated.In Sec.V, we provide a summary and discuss future research directions.

A. Moiré Hamiltonian
Twisted TMD homobilayers with a long-range moiré period has two distinct stacking configurations, of which the twist angle θ between the two layers are, respectively, near 0 • and 180 • .These two configurations are different because each monolayer TMD has a D 3h point-group symmetry without C 2z symmetry (i.e., twofold rotation around out-of-plane ẑ axis).The twisted bilayer with θ close to 180 • can realize a two-orbital Hubbard model on a triangular lattice, as proposed in Ref. 10.
In this work, we focus on valence band states in twisted bilayer WSe 2 with a small twist angle θ near 0 • , motivated by recent experimental studies [14,15].This situation has not been studied much in the theoretical literature so far.As shown in Fig. 1(a), the moiré pattern formed in the twisted bilayer has a period a M ≈ a 0 /|θ|, where a 0 ≈ 3.28 Å is the monolayer lattice constant.In each moiré unit cell (MUC), there are three highsymmetry positions: R M M , R X M and R M X , where M and X, respectively, represent metal and chalcogen atoms, and R β α marks a local position where the α atom in the bottom layer is vertically aligned with the β atom in the top layer.The twisted bilayer has D 3 point-group sym- metry generated by a threefold rotation C 3z around the ẑ axis and a twofold rotation C 2y around the in-plane ŷ axis that swaps the two layers.The D 3 point group is reduced to C 3 when an external out-of-plane displacement field is applied to the system.
In semiconducting TMDs, the topmost valence band states at ±K valleys and Γ valley can be close in energy [16].For small angle twisted bilayer WSe 2 , STM measurement shows that its topmost moiré valence bands originate from ±K valleys instead of Γ valley [14].Therefore, we focus on ±K valleys states.
There is a large valley-dependent spin splitting in the valence bands at ±K valley, which leads to an effective spin-valley locking [11] and reduces the degrees of freedom in the low-energy theory.Therefore, we only consider the spin up (down) valence band in +K(−K) valley, as schematically shown in Fig. 1(c).Furthermore, we treat +K and −K valleys separately in the single-particle Hamiltonian because the two valleys are separated by a large momentum when θ is small [Fig.1(b)].Since the two valleys are related by time-reversal symmetry T , we can focus on +K valley, of which the moiré Hamiltonian is given by [10] where the 2 × 2 matrix is in the layer pseudospin space, the diagonal terms are associated with each layer, and the off-diagonal terms describe the interlayer tunneling.In Eq.( 1), m * is the valence band effective mass, the layer-dependent momentum offset κ ± = [4π/(3a M )](− √ 3/2, ∓1/2) capture the rotation in the momentum space [Fig.1(b)], and ∆ ± (r) is the layerdependent moiré potential given by where V and ψ respectively characterize the amplitude where w is the interlayer tunneling strength.
We take the effective mass m * to be 0.45m 0 following the experimental value [17] of monolayer WSe 2 , where m 0 is the electron rest mass.Other parameters (V, ψ, w) could in principle be estimated using first-principles calculations [9,10,18,19].However, such estimations may suffer from large uncertainties as these parameters are very sensitive to the layer separation that varies spatially in the moiré pattern.Therefore, we treat (V, ψ, w) as phenomenological parameters, and present a systematic study of the moiré band structure as a function of these parameters.At this early stage of the development of the subject, first principles band structure calculations, with their inherent quantitative uncertainties, should be used with caution in developing low energy effective theories with small energy scales, where the relevant band parameters can be obtained from experimental measurements (or can be taken as unknown phenomenological parameters of the effective theory).

B. Layer pseudospin skyrmion
From the continuum Hamiltonian H ↑ in the layer pseudospin space, we can define a scalar potential ∆ 0 and a layer pseudospin magnetic field ∆ as follows We plot the layer pseudospin magnetic field ∆ in Fig. 2(a).The in-plane vector (∆ x , ∆ y ), which accounts for interlayer tunneling, forms vortices and antivortices around R X M and R M X positions, while ∆ z , the z component of ∆, takes maximum and minimum values at these two high-symmetry positions.This spatial profile indicates that ∆ forms a skyrmion lattice, which is characterized by the following winding number N w Here N w is quantized to +1 or −1 depending on the sign of V sin ψ.
In the adiabatic limit where the electron's pseudospin follows the skyrmion texture locally, electron's wave function acquires a real-space Berry phase [20], which can be attributed to an emergent (fictitious) orbital magnetic The effective magnetic flux produced by b z over one MUC is quantized to ±h/e, following Eq.( 5).Fig. 2(b) plots the spatial variation of b z in the moiré pattern, and shows that b z has a strong spatial variation with a large peak value on the order of a few hundreds of teslas, much higher than any real available laboratory magnetic fields.
The skyrmion lattice and the emergent b z field open up the possibility for topological moiré bands.However, we note that the adiabatic limit is not always satisfied in our system, and we find that the skyrmion winding number and the band topology do not have a one-to-one correspondence.
We also define an effective total potential ∆ = ∆ 0 + |∆|.Because the kinetic energy in Eq. ( 1) has a hole type dispersion, low-energy states in our theory are those that are close to the valence band edge.In a semiclassical picture, low-energy states near the band edge tend to be confined near positions where ∆ reaches its maximum value.The maximum positions of ∆ can be at R depending on the exact values of (V, ψ, w), which can have important implications on the band topology, as discussed in the following.

C. Topological phase diagram
We diagonalize the moiré Hamiltonian in Eq. ( 1) using plane-wave expansion based on Bloch's theorem, and show representative moiré band structure in Fig. 3. To discuss band topology, we use C ±K,n to denote the Chern number of the n-th moiré valence band in ±K valleys.
Here we label the moiré valence bands in a descending order of energy, and the topmost moiré valence band in each valley is labelled as the first one.We focus our discussion on +K valley, since C −K,n = −C +K,n because of time reversal symmetry.
We find that the topological character of the moiré bands depends on the precise values of the band parameters.In Fig. 3(a), the first moiré band is topologically trivial with a zero Chern number.By contrast, in Fig. 3(b) with a different set of parameter values, the first moiré band is topologically nontrivial with a finite Chern number.The fact that the topology of the moire bands depends on the details of the parameter values is not surprising since the relevant band Chern number depends on the details of the wavefunction and is not determined uniquely by any symmetry.For the two sets of parameter values used respectively in Figs.3(a) and 3(b), the corresponding skyrmion winding numbers N w are both quantized to +1, which shows that the moiré band topology is not uniquely determined by N w as the adiabatic limit is not always satisfied.
The band topology turns out to have a close connection with the spatial pattern of the effective total potential ∆.For Fig. 3(a), the corresponding ∆ reaches its potential maximum at R M M positions.Therefore, electrons in the first moiré band of Fig. 3 .Therefore, the full system that consists of ±K valleys can realize the Kane-Mele model [22] that includes two time-reversed partner copies of the Haldane model.
To obtain a systematical characterization of the band topology, we present a phase diagram in Fig. 4(a ical even in the parameter space where the first moiré valence band is topologically trivial.This has important experimental consequences since, in principle, these higher topological moiré bands can be studied experimentally if the chemical potential resides in the higher bands.

D. Comparison with STM experiment
To determine which phase has actually been realized in twisted bilayer WSe 2 , we now turn to STM experiments on this system very recently reported in Ref. 14.In this experiment [14], the first LDOS peak at the valence band side (i.e., holes) is found to be primarily localized at R M M positions and the second peak is localized at R X M and R M X positions, which is consistent with the LDOS structure [Fig.6(a)] in the trivial regime of Fig. 4(a).With this comparison between the experiment [14] and our theory, we conclude that the first moiré valence bands in twisted bilayer WSe 2 are most likely topologically trivial.The energy separation ∆E between the first and second LDOS peaks is found to be ∼ 40 meV for twisted bilayer WSe 2 with θ ≈ 3 • in Ref. 14.We plot our theoretical value of ∆E as a function of V and ψ at a fixed value of w in Fig. 6(b).The experimental value ∆E ≈ 40 meV constraints (V, ψ, w) to a finite parameter space that belongs to the topologically trivial regime of Fig. 4(a), but does not lead to a unique determination of (V, ψ, w).We choose a typical set of parameter (V, ψ, w) =(4.4 meV, 5.9, 20 meV), which reproduces the experimental LDOS structure both qualitatively and quantitatively, and use them in all the following calculations.

III. FIELD-TUNABLE LATTICE MODEL
We focus on the first moiré valence band in the topologically trivial regime, and construct an effective tightbinding model for this band in the presence of a layer potential difference V z .Here V z is generated by an external out-of-plane displacement field, and is a tuning knob in controlling the band structure as well as manybody physics.With a finite V z , we replace ∆ ± in the moiré Hamiltonian of Eq. ( 1) by ∆ ± ± V z /2.At V z = 0, the wave function of the first moiré band in +K valley at the two corners of the moiré Brillouin zone, κ + and κ − , are primarily located in the bottom and top layers, respectively.A finite layer potential difference V z shifts the band energies at κ + and κ − in opposite ways, and therefore, can lead to a drastic change in the band structure as demonstrated in Figs.7 (a)-(e).A noticeable effect is that the van Hove saddle points in the band structure can be effectively moved in the moiré Brillouin zone by tuning V z .There is a critical value of V z , at which three van Hove saddle points merge to a single higher-order saddle point [23][24][25] at one of the corners of the moiré Brillouin zone [Fig.7(d)].As a result, the van Hove singularities in the density of states (DOS) can be tuned from below to above half filling by changing V z as shown in Figs.7 (f)-(j).This can have important implications on many-body physics, as discussed in Section IV.
To build up a tight-binding model for this topologically trivial band, we construct localized Wannier states.For this purpose, we choose a gauge such that the bottom layer component of the Bloch wave function at each momentum is real and positive at the origin in real space.A linear superposition of such Bloch states leads to the Wannier state in Fig. 8, which is exponentially localized around the origin (one of the R M M sites) and threefold rotational symmetric.
The corresponding tight-binding model on the triangular lattice formed by R M M sites can be parametrized as where s =↑, ↓ represents spin ↑ and ↓ states associated respectively with +K and −K valleys, R i represents a site in the triangular lattice, and c j,s (c † j,s ) is electron annihilation (creation) operator.t s (R i − R j ) is the hopping parameter, which is constrained by the following relations.(1) Hermiticity of Hamiltonian ( 7) requires that t s (R) = t * s (−R); (2) Three-fold rotational symme- n can also be controlled by V z , as illustrated in Figs.10(c) and 10(d).An important effect is that the phase φ ↑ 1 can be drastically changed by V z .φ ↑ 1 is π at V z = 0, and evolves to 4π/3 (2π/3) when |V z | becomes large enough so that the two layers in the system become effectively decoupled.The dependence of φ ↑ 1 on V z follows the change in the band structure shown in Fig. 7.When the hopping parameters take complex values (i.e., φ ↑ n deviates from 0 or π), they become spin dependent, which leads to effective spin-orbit couplings in the tight-binding model.

IV. HUBBARD MODEL
Many-body interactions are effectively enhanced for electrons in the moiré band with a narrow bandwidth because of the strongly suppressed kinetic energy.By combing the tight-binding Hamiltonian in Eq. ( 7) with electron-electron Coulomb repulsion, we can construct a generalized Hubbard model: where the repulsion U (R i − R j ) between sites i and j is calculated by projecting the Coulomb repulsion Ũ (r) = e 2 /( r) onto the Wannier states.Here is the effective background dielectric constant that can be controlled by the three-dimensional dielectric environment.We take as a free parameter in our theory since its precise value is tunable (and not always precisely known).Numerical values of U 0 (on site repulsion) and U n (n = 1, 2, 3 for repulsion between n-th nearest neighbors) are presented in Fig. 10(b).For a typical value of about 10, the on-site interaction U 0 can be at least one order-of-magnitude greater than the hopping parameters for twist angle θ below 5 • .Therefore, twisted bilayer WSe 2 provides a platform to simulate the generalized Hubbard model on a triangular lattice.Moreover, the hopping parameters can be in situ controlled by an external displacement field.The effective interacting model is a generalized Hubbard model since both interaction and hopping in Eq. ( 8) are not necessarily restricted to being on-site or nearest-neighbor, respectively as the whole many-body Hamiltonian matrix of Eq. ( 8) can be calculated from our moiré band calculations for a given .

A. Heisenberg model
We consider carrier density at half filling, where there is one electron per moiré unit cell in the first moiré valence bands (equivalently, one hole per moiré unit cell when counting from the charge neutrality point of the  9) as a function of φ ↑ 1 .120 • AF ± refer to in-plane antiferromagnetic phases shown in (b) and (c), while FM represents an in-plane ferromagnetic phase.In twisted bilayer WSe2, φ ↑ 1 is constrained between 2π/3 and 4π/3, so only the 120 • AF ± phases are possible.twisted bilayer).The strong on-site repulsion U 0 suppresses double occupation at the same moiré site and gives rise to a Mott insulator.In this Mott limit (where U 0 is very large, much larger than the hopping parameters), the low-energy degrees of freedom are the electron spins at different sites.By retaining only nearestneighbor hopping in the tight-binding model and on-site repulsion U 0 , we can map the Hubbard model in Eq.( 8) to spin Heisenberg model: where the sum over i, j is restricted to nearest neighbors, the prime on the sum indicates that each pair of sites is counted only once, and S i is the spin-1/2 operator at site i.Note that this mapping of the Hubbard model in Eq. ( 8) to the Heisenberg model in Eq. ( 9) involves keeping only the on-site interaction and nearest-neighbor hopping as in the original minimal (rather than the generalized) Hubbard model.In general, a more complete mapping of the full fermion model of Eq. ( 8), i.e., the generalized Hubbard model, to the spin model of Eq. ( 9) is, in principle, possible, but this involves complicated multispin terms beyond the Heisenberg model.This is unnecessary in the current problem since U 0 and t 1 indeed dominates the quantitative physics, thus allowing a mapping from an effective Hubbard model to an effective Heisenberg model of Eq. ( 9).The first two terms in Eq. ( 9) are spin exchange interactions as in a standard anisotropic Heisenberg model, while the last term (S i × S j ) • ẑ describes as an effective Dzyaloshinskii-Moriya (DM) interaction that is generated by the spin-orbit coupling inherent in the tight-binding model of Eq. ( 7).In Eq. ( 9), φ s i,j is the phase of the hopping parameter t s (R i − R j ) between nearest neighbor sites, and ẑ is the unit vector along out-of-plane direction.The relation is used in the simplification that leads to Eq. ( 9).One of the nearest neighbor hopping phases in the spin up channel is shown in Fig. 10 as φ ↑ 1 , and phases for other nearest-neighbor hopping parameters are related to φ ↑ 1 by the three relations given below Eq. (7).Therefore, the single parameter φ ↑ 1 , which is tunable by the layer potential difference V z , determines the ground state of the spin effective Heisenberg model in Eq. (9).
For V z = 0, φ ↑ 1 is π, the DM interaction vanishes, and the model in Eq. ( 9) becomes the standard Heisenberg model with spin SU(2) symmetry on a triangular lattice.This isotropic Heisenberg model with only nearestneighbor exchange coupling has a family of degenerate ground states with three-sublattice 120 • long-range antiferromagnetic (AF) order, which we refer to as 120 • AF states.
For a finite V z , φ ↑ 1 deviates from π, and the finite DM interaction in the spin model ( 9) reduces the spin SU(2) symmetry down to U(1) symmetry, which originates from the valley U(1) symmetry in the Hubbard model of Eq. ( 8).In Fig. 11, we show the calculated classical magnetic phase diagram of Eq. ( 9) as a function of φ ↑ 1 .This diagram is obtained by approximating the spin operator S i as a classical vector with a fixed length and minimizing the energy with Luttinger−Tisza method [26].In our system, φ ↑ 1 takes values between 2π/3 and 4π/3, and crosses π when V z crosses 0, resulting in a sign change in the DM interaction.The DM interaction acts as an anisotropy that favors in-plane spin ordering, and selects a subset of the 120 • AF states to be the ground state.In particular, the spin ground states for V z < 0 and V z > 0 are respectively the 120 • AF − and 120 • AF + phases, which are demonstrated in Figs.11(b) and 11(c).To distinguish these two phases, we choose three sites (A, B, C) along a vertical line in the triangular lattice, as marked in Fig. 11.The spins along the path A → B → C rotate clockwise (anticlockwise) in the 120 • AF − (AF + ) phase.Therefore, these two phases have opposite vector-spin-chirality orders that can be characterized by S A × S B .

B. Mean-field theory
We also perform a mean-field study of the generalized Hubbard model defined by Eq. ( 8) at half filling.The mean-field calculation is not subject to the limit that U 0 |t 1 | and provides an estimation of the charge excitation gap for the interaction-driven correlated insulator at half filling.We use the mean-field ansatz from the spin configuration in the ground state of the Heisenberg model of Eq. ( 9), i.e., the 120 • AF + (AF − ) states for positive (negative) V z .These two different ansatzs for V z > 0 and V z < 0 can also be understood from Fermi surface instability.As shown in Fig. 12, the spin ↑ and ↓ Fermi surfaces in the non-interacting limit have an approximate nesting, with the opposite nesting vector in the momentum space for opposite signed V z .This approximate nesting can lead to interaction-driven instability in the spin density wave (SDW) channel.The SDW order parameter can be taken as c † k+Q,↓ c k,↑ and c † k−Q,↓ c k,↑ , respectively, for V z > 0 and V z < 0.Here we approximate Q by the commensurate wave vector κ + − κ − that connects the two corners of the moiré Brillouin zone.The spin ordering wave vectors ±Q are opposite for opposite signed V z , following the Fermi surface configurations shown in Fig. 12.The order parameters c † k±Q,↓ c k,↑ in momentum space correspond to the 120 • AF ± state in real space.Therefore, the Heisenberg model in the strong coupling limit and the Fermi surface instability in the weak coupling limit are consistent with each other.
With the above mean-field ansatzs, we perform a selfconsistent mean-field calculation for the Hubbard model ( 8) that takes into account hopping up to the third nearest neighbors and the on-site Coulomb repulsion U 0 .In- cluding off-site Coulomb repulsion, which is much smaller than the on-site term U 0 , where i and j are different in the second term of the right hand side of Eq. ( 8) is straightforward, but does not lead to any qualitatively different results, essentially implying a small renormalization of the value of the Hubbard interaction U 0 .The calculated charge gap E G at half filling is shown in Fig. 13.E G is finite for a large range of twist angle θ.Therefore, there is no need to fine tune θ in twisted bilayer WSe 2 in order to realize correlated insulators as is absolutely necessary for twisted bilayer graphene, where the correlated insulator phase is very fragile.We find that E G has a strong dependence on V z particularly for weak interactions (large dielectric constant ).This follows the strong dependence of the non-interacting DOS as well as the nesting degree of Fermi surfaces at half filling on V z , as demonstrate in Fig. 7 and also in Fig. 13(a).A larger non-interacting DOS and a better nesting degree at half filling implies a stronger interaction-driven instability towards symmetry breaking states.As a result, E G can have a dome shape dependence on V z , and the interaction driven insulator at half filling can be turned on and off by V z [Figs.13(b)-(d)].We note that the calculated charge gap in Fig. 13 for reasonable values of are of the order of tens of meVs implying rather robust correlated insulating phases in twisted WSe 2 .

C. Comparison with transport experiment
We compare our theoretical studies with the transport experiment on twisted bilayer WSe 2 in Ref. 15.This experimental paper [15] presents transport study on multiple devices of twisted bilayer WSe 2 with the twist angle θ in the range between 4 • and 5 • , and reports correlated insulators at half filling of the first moiré valence bands.Our theory is consistent with this experiment [15] in key aspects as discussed in the following.
The measured van Hove singularities determined from Hall effect have a strong dependence on displacement field [15].This behavior is captured by our band structure calculation shown in Fig. 7, which shows that the van Hove singularities.
The correlated insulators at half filling develop for a large range of twist angle up to about 5 • , and is controllable via displacement field [15].Our mean-field calculation shown in Fig. 13 provides a qualitative description of this observation.In particular, we also find a dome shape dependence of the insulating gap at half filling on the layer potential difference, as in the experiment [15].This dome-shaped experimental insulating gap is a few (∼2-4) meV typically in Ref. 15 rather than being >10-40 meV or so as we find mostly for the excitation gap in our theory.We note, however, that the measured insulating gap in Ref. 15 is even quantitatively consistent with our theoretical charge gap in Fig. 13 for a large value of (40 or above).This quantitative agreement for large dielectric constant should not be taken too seriously because our mean field theory is bound to overestimate the magnitude of the gap and the experiment measures a transport activation gap which is typically much smaller than the theoretical excitation gap.
The correlated insulating gap at half filling is experimentally found to decrease with increasing out-of-plane magnetic field B z .This is consistent with our in-plane 120 • AF ± states, which turns into canted antiferromagnetic state in the presence of B z .

V. CONCLUSION
In summary, we present a systematic theoretical study of twisted bilayer WSe 2 , and demonstrate the perspective of using this moiré system as a platform to realize interesting single-particle physics as well as many-body physics.For the single-particle moiré bands, we calculate the topological phase diagrams characterized by the valley-contrast Chern numbers.By comparing the theoretical LDOS with STM measurements [14], we conclude that the first moiré valence band is likely to be topologically trivial, whereas the second and third moiré valence bands are likely to be topological.By increasing the hole density in the system, it should be possible to study the topological moiré bands experimentally if one can push the Fermi level into the higher moiré bands.
For the interacting physics, we focus on the first moiré band, and construct a generalized Hubbard model.We show that twisted bilayer WSe 2 can act as a highly tunable Hubbard model simulator.In particular, the layer potential difference V z can drastically change the noninteracting moiré bands, control the charge excitation gap of the correlated insulators at half filling, and generate an effective DM interaction in the corresponding spin Heisenberg model at half-filling.The moiré bands in twisted bilayer WSe 2 are relatively flat over a large range of twist angles θ.Therefore, observation of correlation effects does not require fine tuning of θ in this system, which represents an advantage compared to TBG.We envision that several directions can be explored following our theory.The transport experiment in Ref. [15] has been limited to filling factors within the first moiré valence bands, which are likely to be topologically trivial.It would be interesting to increase the hole doping level, and perform transport study in the second and even third moiré valence bands which likely carry finite valleycontrast Chern numbers.The spin dependent Berry curvatures in these bands can lead to large spin Hall effect.The enhanced Coulomb interactions may drive valley polarization, which, combined with the finite valley Chern number, can lead to quantum anomalous Hall effects.
We show that a field-tunable DM interaction can be realized in the spin Heisenberg model.This DM interaction pins vector spin chirality of the antiferromagnetic ground state.It is desirable to explore effects of DM interaction on spin and magnon transport, and find experimental probes that can distinguish opposite vector spin chiralities.The 120 • AF ± states spontaneously break the U(1) symmetry of the Heisenberg model in Eq. ( 9), which can then support spin superfluidity.
We construct a genearlized Hubbard model on triangular lattice with a field-tunable spin-orbit coupling, and study this model at half filling by mapping it to the Heisenberg model as well as using a mean-field theory.It is conceivable to investigate this model using other techniques and also at other filling factors.Hubbard model on triangular lattice can potentially host a variety of intriguing phases, for example, quantum anomalous Hall insulators [27], chiral superconductors [28], and even spin liquids [29].The inclusion of spin-orbit coupling should enrich the physics.Possible signatures of superconductivity in twisted bilayer WSe 2 have been reported in Ref. [15].Our theoretical model can be a starting point to address exotic many-body physics including superconductivity in this system.

FIG. 1 .
FIG. 1.(a) Moiré superlattices formed in the twisted bilayer.The dots with cyan, red and orange colors indicate, respectively, R M M , R X M and R M X positions with local stacking configurations shown in the insets.The cyan lines mark a moiré unit cell.(b) Brillouin zones associated with the bottom (blue) and top (red) layers, and the moiré Brillouin zone (black).(c) Schematic illustration of band structure in the twisted bilayer.

FIG. 2 .
FIG. 2. (a) The spatial variation of the layer pseudospin magnetic field ∆(r) in the moiré pattern.The arrows represent the x and y components of ∆(r) and the color map shows the z component.The cyan, red and orange dots mark highsymmetry positions as in Fig. 1(a).(b) The effective magnetic field bz(r) that corresponds to the skyrmion field in (a).Parameter values are (θ, V, ψ, w) = (3 • , 5 meV, 0.5π, 20 meV).

FIG. 3 .
FIG. 3. (a) and (b) moiré band structure for different values of (V, ψ).θ is 3 • and w is 20 meV.(c) Local density of states (LDOS) for the moiré bands in (a).The horizontal axis is along a high symmetry line in the moiré pattern.(d) Similar as (c) but for moiré bands in (b).(e) An effective triangular lattice model for the first moiré band in (a).(f) An effective honeycomb lattice model for the first and second moiré bands in (b).In (e) and (f), the cyan, red and orange dots mark R M M , R X M and R M X positions in the moiré pattern.

FIG. 4 .
FIG. 4. (a) Topological phase diagram characterized by the Chern number C+K,1 of the first moiré valence band.(b)The white (green) regime represents parameter space where the potential maximum positions of ∆ are at R M M (R X M /R M X ).The phase boundary between topological (C+K,1 = 0) and trivial (C+K,1 = 0) phases in (a) closely follows the boundary between white and green regimes in (b).

FIG. 5 .
FIG. 5. (a) Topological phase diagram characterized by the Chern number C+K,2 of the second moiré valence band.The colors encode different integer values of C+K,2.(b) Similar as (a) but for the Chern number of the third moiré valence band.

FIG. 6 .
FIG. 6.(a) Local density of states at R M M (blue curve) and R X M (red curve) positions.The first (second) peak marked by blue (red) dashed line is mainly at R M M (R X M ) positions.Parameter values are the same as those for Fig. 3(a).We use ∆E to denote the energy separation between the first and second peaks.(b) ∆E as a function of V and ψ with w = 20 meV and θ = 3 • .The color map shows the value of ∆E in the topologically trivial regime of Fig. 4(a).The white region corresponds to the topological phases of Fig. 4(a), where we do not present the value of ∆E.
(a) are confined to R M M positions, which is verified by the local density of states (LDOS) plotted in Fig. 3(c).It follows that the first band in Fig. 3(a) can be described using a tight-binding model on a triangular lattice formed by R M M sites [Fig.3(e)].As a comparison, the first and second moiré bands in Fig. 3(b) are topological with Chern numbers of −1 and +1, respectively.Electron density in both bands is peaked near R X M and R M X positions [Fig.3(d)], following the potential maximum positions of the corresponding ∆.As shown in Ref. 10, these two topological bands with opposite Chern numbers as a whole can be described by the Haldane model [21] on a honeycomb lattice formed by R X M and R M X sites [Fig.3(f)]

FIG. 7 .
FIG. 7. (a)-(e) The first moiré valence band in +K valley for different values of Vz.The dashed lines mark the Fermi contour at the van Hove energy.Parameter values are (θ, V, ψ, w) = (4 • , 4.4meV, 5.9, 20meV).The plotted band is topologically trivial.(f)-(j) The corresponding density of states for the band shown in (a)-(e).The horizontal axis represents the hole filling factor n/ns.The first moiré valence bands are fully filled (empty) at n/ns =0 (1).
FIG. 8. (a)-(b) Amplitude of W b (r) and Wt(r), which are respectively the bottom and top layer component of the Wannier state.(c)-(d) Phase of W b (r) exp(−iκ+ • r) and Wt(r) exp(−iκ− • r), where the additional phase factors exp(−iκ± • r) make the three-fold rotational symmetry transparent.The black lines mark the effective triangular lattice.Parameter values are the same as those used for Fig. 7(a).

FIG. 10
FIG. 10.(a) |tn| and (b) Un as a function of the twist angle θ. is the effective dielectric constant.Vz is 0 in (a) and (b).(c) |tn| and (d) φ ↑ n as a function of Vz. θ is 4 • in (c) and (d).

FIG. 12 .
FIG. 12. Non-interacting Fermi surfaces at half filling for spin ↑ (blue) and ↓ (red).The dashed vectors indicate an approximate nesting between spin ↑ and ↓ Fermi surfaces.Vz is negative in (a) and positive in (b).

FIG. 13
FIG. 13.(a) The non-interacting density of states at half filling as a function of θ and |Vz|.On the dashed line, van Hove singularities are at half filling.(b) The interaction-driven insulating gap as a function of θ and |Vz|.At a given θ, the gap has a dome shape dependence on |Vz|, which is illustrated in (c) and (d).