Realizing a tunable honeycomb lattice in ABBA-stacked twisted double bilayer WSe 2

,

Theoretically, a strategy to create an ideal emergent honeycomb lattice with sublattice and SU(2) spin symmetries in TMD multilayer structures follows three criteria as first pointed out in Refs.[24][25][26]: (1) The Γ-valley of each TMD layer, where SOC is negligible, appears at the valence band edge (rather than the K/K ′ -valleys), forming an emergent lattice for the doped holes [24][25][26][27].
(2) A symmetry relates the independent high-symmetrystacking MX and XM sites [see Fig. 1(b)], which form the two sublattices of the emergent honeycomb lattice; (3) The energy at the MX and XM sites is higher than the MM site (another high-symmetry stacking location) [28].However, previous papers only proposed these strategies in TMD bilayers.The correct energetics (in criteria (1) and (3)) in these proposals are yet to be confirmed and realized by experiments.For example, Refs.29 and 30 have experimentally identified the Γ-valley physics in twisted bilayer WSe 2 but found it to appear in energies much below the valence band edge, which is located at the K/K ′valleys.
Here, we pursue realizing an ideal emergent honeycomb lattice in ABBA-stacked twisted double bilayer WSe 2 (tdbWSe 2 ), i.e., two AB-and BA-stacked WSe 2 bilayers with a small relative twist angle [see Fig. 1(a)].As experimentally shown in a closely related system, the ABABstacked tdbWSe 2 , the increase in the number of layers enhances the interlayer hybridization driving the Γ-valley to the valence band edge [27,31].However, the ABABstacked tdbWSe 2 does not have the symmetry needed for criterion (2) and only realizes a Γ-valley-based emergent triangular lattice.In this paper, we show that the ABBAstacked tdbWSe 2 can pass the three criteria and realize a tunable ideal honeycomb lattice with both sublattice and SU(2) spin rotation symmetries.
The remainder of this paper is organized as follows.In Sec.II, we provide the microscopic analysis that the ABBA-stacked tdbWSe 2 can pass the three criteria above and, hence, leads to an emergent honeycomb lattice with both sublattice and SU(2) spin rotation symmetries.We construct a continuum model for this tdbWSe 2 system, enabling estimations of the hopping parameters of the emergent honeycomb lattice.In Sec.III, we demonstrate the intriguing tunability of this system by studying orders induced by an in-plane magnetic field at the hole filling ν = 2, namely, the half-filling of the emergent honeycomb lattice.The phase diagram is obtained using the Hartree-Fock analysis of the extended Hubbard model on the emergent honeycomb lattice.Finally, in Sec.IV, we summarize our results and experimental signatures of the predicted orders.

II. MODELING THE EMERGENT HONEYCOMB LATTICE
The ABBA-stacked tdbWSe 2 can satisfy the three criteria above for the following reasons.First, this structure's increased number of layers leads to stronger interlayer tunnelings compared to a single bilayer, pushing the Γ-valley to the top of the valence bands [32,33].Second, inherited from the D 3h point group symmetry of each WSe 2 layer, the ABBA-stacked tdbWSe 2 manifests a three-fold rotation C 3z around the z-axis, and a twofold rotation C 2y around the in-plane y-axis, where the y-axis is chosen to align with the xy-plane projection of a bond between the tungsten atom and the selenium atom in the monolayer WSe 2 .These crystal symmetries ensure a sublattice-symmetric honeycomb lattice composed of the MX and XM sites.Third, from the ab initio tightbinding model in Appendix A, we find the correct energy hierarchy amongst the MX, XM, and MM sites to ensure an emergent honeycomb lattice.
To estimate the hopping parameters in the emergent honeycomb lattice based on the three criteria, we need to obtain the dispersion of the two topmost moiré valence bands.Adapting the approach in Ref. [25,31], we begin with a four-layer continuum Hamiltonian for our twisted ABBA-stacked system that preserves C 2y , C 3z , time-reversal and SU(2) spin rotation symmetries (due to negligible SOC near the Γ-valley [34]): (1) H 4L is identical for both spin species.
The intralayer potentials take the forms ∆ 1 (r) = ∆ 4 (r) = V 1 for the first and fourth layer, and ∆ 2 (r) =  We now estimate the parameters in the intralayer potentials and interlayer tunnelings in Eq. ( 1).These potentials and tunnelings in the tdbWSe 2 evaluated at the MM, and MX/XM sites should be fitted according to the band structures of the untwisted ABBA-stacked double bilayer WSe 2 in the three high-symmetry stacking configurations MM, and MX/XM [see Appendix A for more details] [36].The angle-independent parameters for the potentials and tunnelings are estimated as follows: 2 ) = (200, −159, −8) meV, ϕ = −0.17,and (V 12 , V 23 ) = (184, 356, −9) meV.We comment that these parameter estimates are obtained based on assumptions on the interlayer distances in the double bilayer WSe 2 detailed in App. A. The assumed interlayer distances produce a representative example of tdbWSe 2 satisfying the energetics criteria allowing us to estimate bandwidth and demonstrate the tunability of the emergent honeycomb lattice, which are the main purposes of our work.With the above parameters, we can solve Hamiltonian (1) for the moiré valence bands of tdbWSe 2 .For example, Figure 1(d) shows the band structure along the γ-m-κ-γ path in the moiré Brillouin zone (mBZ) at a twist angle of θ = 2 • .From the two topmost moiré valence bands, especially the presence of a Dirac cone between them, we confirm that ABBA-stacked tdbWSe 2 with a hole filling less than 4 per moiré unit cell captures an emergent honeycomb lattice with both sublattice and SU(2) spin rotation symmetries.At a twist angle of θ = 2 • , the bandwidth is around 10 meV (which is much smaller than the ∼2 eV bandwidth of graphene [37]).
At hole filling ν = 2, i.e., two holes per moiré unit cell, the Fermi level is precisely at the Dirac point protected by C 3z and C 2y symmetries.This results in a semimetal (non-interacting) ground state, as shown in Fig. 1(d).The narrow bandwidth of the tdbWSe 2 -realized honeycomb lattice allows for an intriguing way to tune the electronic structure: a realistic in-plane magnetic field can "dope" the semimetal and create finite Fermi surfaces located around both κ and κ ′ points of the mBZ for both spin species.Note that such tuning is impossible for the K/K ′ -valley-based TMD moiré system due to the strong SOC therein.The "doping" under the inplane magnetic field is the result of the Zeeman effect.The thinness of the double-bilayer structure allows us to ignore the orbital effect of the in-plane field.Without loss of generality, we assume the in-plane magnetic field B is applied along the x direction, resulting in a Zeeman energy of h x = 1 2 gµ B B, where µ B is the Bohr magneton (5.79 × 10 −2 meV/T,) and the g-factor is roughly 2 due to the negligible SOC near the Γ-valley.For example, at a twist angle of θ = 2 • , a magnetic field of 13 T, typically The schematic picture of the ABBA-stacked tdbWSe2.Within each of the AB-and BA-stacked bilayer WSe2, the orientation of the two layers differs by 180 • .The two bilayers are twisted so that the middle two layers have a small twist angle θ between them.(b) The real-space moiré pattern created by the twisted central two layers (blue and orange) with the darker (lighter) for the metal (chalcogenide) atom.The three high-symmetry-stacking sites in one moiré unit cell (black dashed lines) are MM (brown, metal atoms on one layer aligned with metal atoms on the adjacent layer), XM (blue, chalcogenide atom on one layer aligned with metal atoms on the adjacent layer), and MX (grey, metal atom on one layer aligned with chalcogenide atoms on the adjacent layer).Here, the alignment of atoms refers to the central two layers in the double bilayer structure; (c) Moiré Brillouin zone with four high-symmetry points γ, m, κ, and κ ′ .G1,.. achievable experimentally, creates a Zeeman splitting energy of 0.8 meV.This Zeeman splitting leads to a finite Fermi surface of each spin species roughly halfway in energy between the Dirac point (at the κ point in the mBZ) and the van Hove singularity (at the m point), as shown in Fig. 2. In the presence of the in-plane field, a magnetic particle-hole instability that opens a charge gap is expected even at weak interactions because of the nesting between the Fermi surfaces of the two spin species.This instability would lead to a magnetic order that spontaneously breaks the remaining U(1) spin rotation symmetry about the x-axis, namely, the in-plane field direction.
Conceptually, a similar in-plane-field-induced particlehole instability applies to graphene.The previous work [3] shows that the resulting state is a canted Néel antiferromagnetic (AF) insulator, whose magnetizations pro-jected onto the plane perpendicular to the applied field are opposite in the two sublattices of the honeycomb.However, realizing such a state in graphene requires a very high in-plane field of 10 2 ∼ 10 3 T [3].As we show below, a similar canted Néel AF insulator can be found in the tdbWSe 2 -realized honeycomb lattice, whose bandwidth is engineered around the meV scale, at a realistic in-plane field of a few Tesla.

III. EXTENDED HUBBARD MODEL AND PHASE DIAGRAM UNDER IN-PLANE FIELD
Aside from the bandwidth distinction with graphene, the band structure of the tdbWSe 2 -realized honeycomb lattice exhibits a more pronounced contribution from the next-nearest-neighbor (NNN) hopping t 2 on the honeycomb lattice.The effect of t 2 is manifested by the asymmetry (about the Dirac cone) between the two topmost valence bands shown in Fig. 1 (d).In contrast to the graphene-inspired model of [3] with only the nearestneighbor (NN) hopping t 1 and the on-site Hubbard interaction, we consider an extended Hubbard model with longer-range hoppings and interactions pertaining to our tdbWSe 2 -realized honeycomb lattice system.Our Hamiltonian consists of three terms: the hopping H h , the Zeeman energy H Z from the in-plane magnetic field, and the extended Hubbard interaction H int .The term Here, c iσ is the fermion operator of the hole with spin σ at the site i of the emergent honeycomb lattice.Both t 1 and t 2 are real because of the SU(2) spin rotation symmetry and the time-reversal symmetry.t 1,2 are are tunable by the twist angle θ.For instance, at θ = 2 • , we estimate t 1 = 1.8 meV and t 2 = 0.2 meV based on the moiré band structure obtained from the continuum model Eq.( 1).The ratio t 2 /t 1 varies as we change the twist angle: 0.018 for θ = 1 • , 0.11 for θ = 2 • , and 0.12 for θ = 3 • .In the search for possible phases below, we treat t 2 /t 1 as a parameter to explore the effect of longer-range hopping.The Zeeman energy H Z is given by where the spin quantization axis conveniently is chosen along the x direction.Under H Z , the spin symmetry is reduced to a U(1) spin rotation in the yz plane.We fix h x = 0.4t 1 in the following.For the twist angle θ = 2 • , h x = 0.4t 1 (corresponding to a 13 T in-plane magnetic field) sets the Fermi level roughly halfway between the Dirac point at the κ/κ ′ point and the van Hove singularity at the m point of the mBZ.
The range of interactions in TMD-based superlattices typically extends beyond on-site interactions.Here, we consider the extended Hubbard interaction with the on-site interaction strength U 0 , and the NN Hubbard interaction of strength U 1 .In the following study of the phase diagram, we treat U 0 and U 1 as two free parameters to explore their effect on the canted Néel AF insulator and other phases.
We start with understanding the effect of the NNN hopping t 2 , and present the ground-state phase diagram obtained using the Hartree-Fock method in Fig. 3(b).In this phase diagram, we consider a parameter range U 0 ∈ [t 1 , 9t 1 ] and t 2 ∈ [0, 0.6t 1 ] while setting U 1 = 0. Using the numerical Hartree-Fock method, we find three possible phases: the canted Néel AF insulator, canted √ 3 × √ 3 spin density wave (SDW), and a metal.At hole filling ν = 2, the Fermi surfaces of the two spin species are always perfectly nested for t 2 ≲ t 1 /3 at both the κ and κ ′ points of the mBZ.This perfect nesting results in a Fermi surface instability towards a robust canted Néel AF insulator, spontaneously breaking the U(1) spin rotation symmetry in the yz plane as mentioned above.Beyond t 2 > t 1 /3, the system's ground state is metallic at small U 0 .The metallic state is due to an extra Fermi surface of the spin-down hole appearing around the γ point when t 2 > t 1 /3.Consequently, the Fermi surfaces of the two spin species are no longer nested at ν = 2.As U 0 increases, the system can enter an insulating SDW phase with an enlarged √ 3 × √ 3 unit cell.This SDW phase exhibits a 120 • AF order in the yz plane on each sublattice of the honeycomb.The chiralities of the two 120 • AF orders are opposite.For the tdbWSe 2 -realized honeycomb lattice, t 2 /t 1 < 1/3 for twist angle θ < 4 • .Hence, restricted to only on-site interactions, the canted Néel AF insulator is the most relevant in-plane-field-induced phase for this material at hole filling ν = 2.
In addition to the NNN hopping t 2 , we also study the effect of the NN Hubbard interaction U 1 .We present the phase diagram obtained from the Hartree-Fock calculation as a function of U 0 and U 1 as shown in Fig. 3(c).For this phase diagram, we set t 2 = 0 for simplicity.The canted Néel AF insulator is robust at small U 1 .As U 1 increases, the system is driven to an insulating phase with the charge distribution polarized to one of the sublattices, spontaneously breaking the system's C 2y symmetry.In contrast to the canted Néel AF insulator, the yz-plane U(1) spin rotation symmetry is preserved in this phase.The phase boundary between the AF insulator and the sublattice-polarized insulator is roughly at x .This is the exact phase boundary in the strong interaction limit (where the hopping terms H h are neglected).In the tdbWSe 2 -realized honeycomb lattice, the screening of the Coulomb interaction determines the relative strength between U 0 and U 1 .Stronger screening favors the canted Néel AF insulator, while weaker screening favors the sublattice-polarized insulator.In Appendix B, we present an analytical mean-field analysis in the weak-coupling limit U 0,1 ≪ t 1 and find that the canted √ 3× √ 3 SDW insulator can also compete with the canted Néel AF insulator when U 1 /U 0 increases.

IV. SUMMARY AND EXPERIMENTAL IMPLICATION
In summary, we show that the ABBA-stacked tdbWSe 2 can serve as a realistic and tunable platform to simulate Γ-valley honeycomb lattice with both sublattice and SU(2) spin rotation symmetries.We develop a continuum model to estimate the relevant hopping parameters of this honeycomb lattice model, including the NN and the NNN hoppings t 1,2 .The small bandwidth, found to be at the meV scale, and the ratio t 2 /t 1 are both tunable via the twist angle θ.We show that the small bandwidth enables the effective control of the electronic structure via an in-plane magnetic field.At a hole filling ν = 2, an in-plane field of a few Teslas can dope this system from a Dirac semimetal to a metal having finite Fermi surfaces with instabilities.To demonstrate this tunability and understand the resulting instabilities, we construct an extended Hubbard model for this honeycomb lattice and perform a numerical Hartree-Fock calculation for the ground-state phase diagram at ν = 2.We find a robust canted Néel AF insulating phase when the on-site Hubbard interaction U 0 dominates.This phase is conceptually similar to the magnetic-field-induced canted AF insulator studied in the context of graphene by Ref. [3].However, our finding resides in a more realistic parameter regime.Moreover, a competing sublatticepolarized insulator is found in our extended Hubbard model when the NN Hubbard interaction U 1 increases beyond a threshold.
The results presented in this work have direct experimental implications.First, for ABBA-stacked tdbWSe 2 at hole filling ν = 2, an in-plane magnetic field is predicted to induce an insulating ground state.There are two possible insulating ground states: the canted Néel AF insulator and the sublattice-polarized insulator.The former can be detected through its spin texture perpendicular to the in-plane field direction, potentially using spin-polarized scanning tunneling microscopy [38].The finite-temperature phase transition into the canted Néel AF order is a Berezinskii-Kosterlitz-Thouless transition associated with the remaining U(1) spin rotation symmetry under the in-plane field.This metal-insulator transition can be accessed by lowering the temperature at a given in-plane field or by increasing the in-plane field at a low temperature starting from the Dirac semi-metal.To be more precise, at zero temperature (and in the absence of disorder), there should be no threshold for the interaction strength or in-plane magnetic field -as long as they are both finite, the insulator phase exists.However, at a finite temperature, one has to turn on an in-plane field and interaction strength beyond certain thresholds to observe the insulating phase with the spin orders.The thresholds can be estimated by comparing the temperature with the gap sizes of the insulating phase at zero temperature.Therefore, at sufficiently low-temperature, one can access a metal-insulator transition using an inplane magnetic field.The sublattice-polarized insulator has no non-trivial spin texture (other than the magnetization along the in-plane field direction).Instead, the spontaneous charge polarization in the two sublattices of the emergent honeycomb lattice leads to a finite electric polarization in the tdbWSe 2 sample along the z direction.Consequently, a signature of this sublatticepolarized insulating phase is the hysteresis of the z-axis electric polarization as a function of the z-direcction electric field.The finite-temperature phase transition into the sublattice-polarized phase is in the Ising universality class.Switching between the two insulating phases can be implemented by tuning the ratio between U 0 and U 1 , which is achievable via adjusting the twist angle and/or sample's distance from the gates.
To show that criteria (1) and ( 3) introduced in the main text can be satisfied by the tdbWSe 2 , the first task is to obtain its band structure.However, it is not easy to derive the band structure of tdbWSe 2 directly at the atomic level because the emergent moiré superlattice contains up to ten thousand atoms even within just one moiré unit cell.One way to approximate the band structure of the twisted system is to use a continuum model with parameters fitted according to the untwisted system in the high-symmetry stacking configurations (MM, MX, and XM) that appear in the moiré unit cell of the twisted system [36].We will provide the details later in this section.Therefore, in this section, we will first discuss the method to obtain the band structure of untwisted systems.
We calculate the band structure of the untwisted system using the ab initio tight-binding Hamiltonian [34].The logic is to first construct a tight-binding model for the monolayer by considering five d orbitals of one tungsten atom and six p orbitals of the two selenium atoms in the monolayer unit cell, and including spinorbit coupling over the entire monolayer Brillouin zone (BZ) (i.e., not just the small mBZ).It is called ab initio because hopping parameters are evaluated from the density functional theory.After constructing the tightbinding model for the monolayer, we can model the interlayer tunneling between the adjacent layers by considering the hopping between the two selenium atoms on the adjacent layers.Combining these considerations, we can construct the Hamiltonian for the (untwisted) double bilayer WSe 2 .Finally, we diagonalize the matrix of size 88 by 88 (11 orbitals × 2 spins × 4 layers) at each momentum to obtain the band structure E r (k), where r = {MM, XM, MX} represents the highsymmetry stacking configuration, and k is the momentum in the entire monolayer BZ.For the full recipe and the numerical details, we refer to Ref. 34.
Our goal is to show that tdbWSe 2 can satisfy criteria (1) and (3) so that a honeycomb lattice with both sublattice and SU(2) spin rotation symmetries emerges at the moiré scale.That is to say, we require E r (Γ) > E r (K) for r =MM and MX/XM, and E MX (Γ) = E XM (Γ) > E MM (Γ).Obtaining these energies from the model above requires us to make a choice for the interlayer distances in the tdbWSe 2 structure, which are unknown parameters.Our main goal is to demonstrate the possibility and the tunability of the emergent honeycomb lattice from tdbWSe 2 and estimate the relevant energy scale.We will not attempt to calculate the relaxed interlayer distances.Instead, we will make assumptions on them based on the following reasoning.In light of the C 2y symmetry which flips the layers, the distance d 12 between the first and second layer should be the same as the distance d 34 between the third and fourth layer.For the distance between the central two layers d 23 , because the MM stacking corresponds to the case where all tungsten atoms (also for selenium atoms) on both layers are stacked with perfect alignment, respectively, the repulsion from the interlayer interaction should be the strongest.Therefore, the value of d 23 in the MM stacking configuration should be larger than d 12 and d 23 in the MX/XM stacking configuration.With these considerations, we perturb the interlayer distances d 12 , d 23 for the MM stacking, and d 23 for the MX/XM stacking, around the interlayer distance from the bulk WSe 2 data, ∼ 6.5 Å [35], and find that the system can satisfy the energetics hierarchy in criteria ( 1) and ( 3) in a large parameter space around d 12 =6 Å.As an example that passes criteria (1) and (3), we choose the following interlayer distances (d 12 , d 23 , d 34 )=(6,6.25,6)Å for MM stacking and (d 12 , d 23 , d 34 )=(6,5.8,6)Å for MX and XM stacking for the estimation of the parameters in the continuum model.Note that this choice has smaller layer distances than the bulk data.In principle, one can consider adding pressure to reduce the interlayer distance if needed.For example, a range of d 23 from 6 to 6.5 Å (from 5.7 to 6.3 Å) for MM (MX and XM) stacking with d 23 = 5.8 Å (d 23 = 6.25 Å) for MX and XM (MM) stacking fixed can also satisfy the both criteria.Now, we present the valance band structure for the untwisted MM stacking in Fig. 4(a), and untwisted MX/XM stacking in Fig. 4(b).The corresponding valence band edge energies at the Γ and K valley for these two stackings are listed in Table I, which validate the two criteria.
To model the twisted moiré system from the untwisted system, we can represent the energy at the highsymmetry-stacking sites in one moiré unit cell using the energies in an untwisted system with the same type of stacking configuration.Then we interpolate them using up to the first harmonics of the moiré periodicity [36].Therefore, by extracting the 4 topmost valence bands at the Γ-valley for the untwisted systems in the MM and MX/XM stacking configurations, we can estimate the unknown parameters in the moiré continuum Hamiltonian Eq. ( 1) in the twisted system.To do that, we first will drop the kinetic term in Eq. ( 1), as it is near Γ valley, and then diagonalize the potential term with r =MM and r =MX/XM, respectively, which gives four discrete energies for each site.The fitting procedure is to adjust these four energies at r =MM or r =MX/XM to match the energies obtained from the four valence bands in the ab initio tight-binding Hamiltonian at Γ valley, i.e., the first two columns in Table I.In this appendix, we present the analytical mean-field study of the Fermi surface instabilities in our extended Hubbard model (with the in-plane magnetic field) in the weak-coupling limit U 0,1 ≪ t 1 at hole filling ν = 2.We find that, at zero temperature, dominant instabilities include both the canted Néel AF order and the canted √ 3 × √ 3 SDW order.Their competition is controlled by the ratio U 1 /U 0 .We also discuss the relationship between our weak-coupling analysis and the Hartree-Fock phase diagrams in Fig. 3 obtained in the stronger coupling regime U 0 ≳ t 1 .
We start with the noninteracting part of the Hamiltonian H 0 = H h + H Z containing the hopping term and the Zeeman energy.At hole filling ν = 2, for small t 2 , h x compared to t 1 , the Fermi surfaces are centered at the κ and κ ′ points of the mBZ.We expand the dispersion around the κ/κ ′ point: Here, k is the momentum measured from the κ/κ ′ point.The dispersion is expanded to the linear order of k, which suffices to study the Fermi surface instabilities.The Fermi velocity v ⊺ is the fermionic spinor in the sublattice space at valley τ = κ/κ ′ with spin σ =↑, ↓ along the in-plane field direction x. ⃗ Λ comprises the Pauli-X and Y matrices acting on the sublattice space.Note that t 2 's contribution only starts to appear in the order-k 2 terms of the dispersion.This contribution is independent of valley, spin, and sublattice.Hence, a small t 2 is unimportant for the Fermi surface instability in the weak-coupling limit.
At non-zero h x , the Fermi surfaces of the spin-up fermions are perfectly nested with those of the spindown fermions.Consequently, particle-hole instabilities towards magnetic ground states with spontaneous hybridization between the two spin species are expected.The magnetic orders in these ground states spontaneously break the U(1) spin rotation symmetry in the yz plane (perpendicular to the in-plane field direction).In the limit h x → 0, the instabilities towards such magnetically-ordered ground states vanish along with the Fermi surfaces in the noninteracting band structure.Hence, these magnetic orders, to be discussed in detail below, are dubbed in-plane-field-induced magnetic orders.
At hole filling ν = 2 (and with finite h x > 0), we begin our analysis of possible instabilities by writing down the noninteracting bands crossing the Fermi level: for the two valleys τ = κ, κ ′ and the two spin species σ =↑, ↓.The energies E ±,σ (k) are given by E ±,σ (k) = h x σ ±v F |k| with ± playing the role of a band index.This band index ± is locked with the spin index σ for the bands that cross the Fermi level.The fermion operators c τ,±,σ are defined as and where θ(k) is the azimuth angle of k.The spin-up (spindown) fermions have two particle-like (hole-like) Fermi surface pockets, one around the κ point and the other around the κ ′ point.Hence, possible particle-hole instabilities (at weak couplings) can occur through either intervalley or intravalley hybridizations of the two spin species.
Here, τ represents the opposite valley of τ .The mean-field decomposition of the extended Hubbard interaction corresponding to the intervalley hybridization reads Here, σ represents the opposite spin of σ.N is the system size.The band index is suppressed in the shorthands c τ,↑ (k) = c τ,−,↑ (k) and c τ,↓ (k) = c τ,+,↓ (k) due to the locking of the spin and the band indices for the bands that cross the Fermi level.In the derivation of H inter int , the terms irrelevant to the intervalley hybridization of opposite spin species and the terms of linear or higher order in k are dropped.Note that the first (second) term in the quadratic terms in the first line of Eq. (B5) only couples to the p-wave (s-wave) channel of the order parameter c † τ,−,↑ (k)c τ ,+,↓ (k) .We can consider the s-wave channel (with orbital angular momentum |l| = 0) and the p-wave channel (with orbital angular momentum |l| = 1) separately.At the mean-field level, the saddle point equations in the s-wave channel only depend on U 1 , while those in the p-wave channel only depend on U 0 .
For the intravalley hybridization associated with the mean-field order parameter c † τ,−,↑ (k)c τ,+,↓ (k) , the mean-field decomposition of interaction terms reads Similar to the intervalley case, based on the form of H intra int , we only need to consider the s-wave and the p-wave channel of the order parameter c † τ,−,↑ (k)c τ,+,↓ (k) .Also, it is natural to expect that the order parameters in the two valleys τ = κ, κ ′ share the same orbital angular momentum.With this expectation, we find that the saddle point equations in the s-wave channel only depend on U 0 while those in the pwave channel only depend on U 1 .
At zero temperature, we find the two most dominant mean-field solutions with the largest gaps opened around the Fermi surface.The corresponding orders are then expected to have the highest mean-field transition temperatures.One of these mean-field solutions is given by the intervalley hybridization in the s-wave channel with c † κ,−,↑ c κ ′ ,+,↓ = c † κ ′ ,−,↑ c κ,+,↓ , which gives rise to the canted √ 3× √ 3 SDW order, as illustrated in Fig. 3(b).In this phase, the gap opened at the Fermi surface is given by 2h x exp − 2 3U1ρ , where the density of states ρ can be approximated to a constant 4hx 3πt 2 1 within the energy interval between two Dirac points (at energies −h x and h x relative to the Fermi level).The other dominant mean-field solution is given by the intravalley hybridization in the i • r − ϕ) for the second and third layer (the opposite sign of ϕ as a result of C 2y symmetry).The interlayer tunnelings take the forms ∆ 12 (r) = ∆ 34 (r) = V 12 and ∆ 23 (r) = V cos(G i • r).Here, FIG. 1.(a)The schematic picture of the ABBA-stacked tdbWSe2.Within each of the AB-and BA-stacked bilayer WSe2, the orientation of the two layers differs by 180 • .The two bilayers are twisted so that the middle two layers have a small twist angle θ between them.(b) The real-space moiré pattern created by the twisted central two layers (blue and orange) with the darker (lighter) for the metal (chalcogenide) atom.The three high-symmetry-stacking sites in one moiré unit cell (black dashed lines) are MM (brown, metal atoms on one layer aligned with metal atoms on the adjacent layer), XM (blue, chalcogenide atom on one layer aligned with metal atoms on the adjacent layer), and MX (grey, metal atom on one layer aligned with chalcogenide atoms on the adjacent layer).Here, the alignment of atoms refers to the central two layers in the double bilayer structure; (c) Moiré Brillouin zone with four high-symmetry points γ, m, κ, and κ ′ .G1,..,6 are the moiré reciprocal vectors; (d) Moiré band structure obtained from the continuum model at twist angle θ = 2 • along the path γ-m-κ-γ [blue dashed line shown in (c)].The red dashed line indicates Fermi energy at hole filling ν = 2 (two holes per morié unit cell).

FIG. 2 .
FIG. 2.The two topmost moiré valence bands shown in Fig.1(d) split by the Zeeman energy shift due to an in-plane magnetic field of 13 T.The arrows indicate the spin along the in-plane field direction.The red dashed line indicates the Fermi energy at hole filling ν = 2, i.e., two holes per moiré unit cell and half-filling for the emergent honeycomb lattice.
FIG. 3. (a) The color wheel shows the six spin orientations in the yz plane.(b) The phase diagram as a function of U0 and t2.(c) The phase diagram as a function of U0 and U1.Here, the size of the dot indicates the number density of each site, while the color indicates the spin orientation in the yz plane as shown in the color wheel.The black dots in the sublattice polarized state indicate zero spin polarization.

TABLE I .
First four valence band (from the top to bottom) energies (eV) at Γ and K-valley for MM and MX/XM stackings