{'fraction_non_alphanumeric': 0.03462311026535574, 'fraction_numerical': 0.01321492758219685, 'mean_word_length': 4.5742486741308195, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Realization of Nearly Dispersionless Bands with Strong Orbital Anisotropy from Destructive Interference in Twisted Bilayer MoS2

Recently, the twist angle between adjacent sheets of stacked van der Waals materials emerged as a new knob to engineer correlated states of matter in two-dimensional heterostructures in a controlled manner, giving rise to emergent phenomena such as superconductivity or correlated insulating states. Here,we use an ab initio based approach to characterize the electronic properties of twisted bilayer MoS2. We report that, in marked contrast to twisted bilayer graphene, slightly hole-doped MoS2 realizes a strongly asymmetric px-py Hubbard model on the honeycomb lattice, with two almost entirely dispersionless bands emerging due to destructive interference. We study the collective behavior of twisted bilayer MoS2 in the presence of interactions, and characterize an array of different magnetic and orbitally-ordered correlated phases,which may be susceptible to quantum fluctuations giving rise to exotic, purely quantum, states of matter.


Introduction
Two-dimensional van der Waals materials constitute a versatile platform to realize quantum states by design, as they can be synthesized in many different stacking conditions (1), offer a wide variety of chemical compositions and are easily manipulated by back gates, strain and the like. Stacking two sheets of van der Waals materials atop each other at a relative twist has recently emerged as a vibrant research direction to enhance the role of electronic interactions, with first reports on twisted bilayer graphene (2)(3)(4)(5)(6) and other van der Waals materials stacked atop each other at a twist (7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17) displaying features of correlated physics that afford an unprecedented level of control. In particular, bi-, tri-and quadruple-layer graphene as well as twisted few-layer transition metal dichalcogenides (TMDs) are currently under intense experimental scrutiny (13,(18)(19)(20)(21)(22)(23)(24)(25)(26). By forming a Moiré supercell at small twist angles, a large unit cell in real space emerges for twisted systems, which due to quantum interference effects leads to a quasi-two-dimensional system with strongly quenched kinetic energy scales. This reduction in kinetic energy scale, signaled by the emergence of flat electron bands, in turn enhances the role of electronic interactions in these systems.
Whereas the flatting of band dispersions in two-dimensional Moiré superlattices results mainly from the localization of charge density distributions by the Moiré potential, a wellknown alternate pathway to flat bands can occur in certain lattices such as the Lieb and the Kagome lattices, where geometric considerations permit the formation of perfectly-localized electronic states on plaquettes and hexagons, respectively, that are eigenstates of the kinetic Hamiltonian due to destructive interference between lattice hopping matrix elements (27). Such flat band systems can give rise to many interesting phenomena, such as the formation of nontrivial topology when time reversal symmetry is broken, or other exotic quantum phases of matter due to their susceptibility to quantum fluctuations and electronic correlations (28).
Here, we demonstrate that both flat-band mechanisms can be engineered to coexist in twisted bilayers of MoS 2 (tbMoS 2 ); a TMD of direct experimental relevance that has been extensively studied from synthesis to applications (29,30). We confirm that families of flat bands emerge when two sheets of MoS 2 are stacked at a twist (12,31) due to Moiré potentials. Our large-scale ab initio based simulations show that while the first set of engineered flat bands closest to the edge of the band gap with twist angles close to Θ ≈ 0 • can be used to effectively engineer a non-degenerate electronic flat band in analogy to a single layer of graphene at meV energy scales, more intriguingly, the next set of flat bands instead realizes a strongly asymmetric flat band p x -p y honeycomb lattice (32,33). Both of these families of bands should be accessible experimentally via gating. The strongly asymmetric nature of this p x -p y honeycomb lattice is in marked contrast to the much discussed case of twisted bilayer graphene, where an approximately symmetric version of such a Hamiltonian is now believed to describe the low-energy flat band structures found at small twist angle (34)(35)(36)(37)(38). The strongly asymmetric p x -p y honeycomb model realized here features two almost entirely dispersionless flat bands that touch the top and the bottom of graphene-like Dirac bands at the Gamma point, respectively. These flat bands originate from destructive interference, in analogy to flat bands in the Lieb and the Kagome lattices (27), and will be referred to as ultra-flat bands in the following discussion, to distinguish them from other bands with quenched kinetic energy scales. In addition, these ultra-flat bands are topologically non-trivial once time-reversal symmetry is broken via spin-orbital cou-pling (39). Previously, the p x -p y model was studied in the context of cold gases where exotic correlated phases were predicted (32,40,41), as well as in semiconductor microcavities (42) and certain 2D organometallic frameworks (43,44). Our findings elevate tbMoS 2 to a novel platform where effects of ultra-flat bands can be studied systematically in a strongly-correlated solid-state setting. Notably, in the strong-coupling regime, the p x -p y model amended by Hubbard and Hund's interactions gives rise to a spin-orbital honeycomb model which -depending on the specific parameters and symmetries of the model -hosts magnetic, orbital as well as valence-bond orderings, or even more exotic quantum spin-orbital liquid phases (45)(46)(47).
With this, our work adds an unprecedented type of lattice model -the highly asymmetric p x -p y Hubbard model -to the growing list of systems that can effectively be engineered using the twist angle between multiple layers. This is particularly intriguing as we maintain the full advantages that come with two-dimensional van der Waals materials, such as relative simplicity of the chemical composition and controllability of the material properties; e.g. of the filling (by a back gate), electric tunability (by displacement fields) or the band width of the model (by the twist angle).

ab initio characterization of twisted MoS 2
We first characterize the low-energy electronic properties of twisted bilayer MoS 2 using density functional theory (DFT) calculations (see Methods). DFT in particular has established itself as a reliable tool to provide theoretical guidance and to predict the band structure of many twisted bi-and multilayer materials (8,13,15). However, such a first principles characterization becomes numerically very demanding as the twist angle Θ approaches small values and the unit cell becomes very large entailing many atoms (of the order of a few thousands and more). Nevertheless, it is that limit in which strong band-narrowing effects and as a consequence prominent  we show the ab initio band structure of the twisted material after relaxation, where we find two families of bands which will become increasingly flat and start to detach from all other bands, as the twist angle is lowered. We mark these bands by blue and red color in panel (b), which shows results for decreasing angles from Θ = 3.16 • to 2.28 • . The bandwidth of these two energetically-separated groups of bands is summarized in panel (c) of Fig. 1. We find that the bandwidth of these two bands shrinks drastically as the angle is decreased, yielding band widths of the order of 10 meV as the angle approaches Θ ≈ 2 • . Note that these flat bands near the top of the valence bands originate from the states around the Γ point in the Brillouin zone of the primitive unit cell of untwisted MoS 2 , with both S p z and Mo d z 2 characters. This is different from the case of twisted WSe 2 , in which the top valence flat bands originate from the states around the K point in the Brillouin zone of the primitive unit cell (13). Since in other TMDs, such as MoSe 2 and WS 2 , the top of the valence band in the untwisted bilayer is also located at the Γ point in the Brillouin zone (48,49), the physics we discussed here transfers also to those materials being twisted.
The upper bands in Fig. 1 (marked in blue) show a Dirac cone at the K point and behave very similar to the bands found for monolayer graphene (with the exception of a reduced band width).
They are spin degenerate in nature, but feature no additional degeneracy except at certain high symmetry points. Instead, the next set of bands (marked in red) is essential to our work. They too feature a Dirac cone at the K point, but also feature two additional ultra-flat bands at the top and bottom in addition to a band structure similar to graphene. The ratio between the width of the ultra-flat and the flat bands decreases as the angle is decreased, but saturates in our calculations as a twist-angle of Θ ≈ 2.28 • is approached. We attribute this saturation to lattice relaxation effects; note however that the overall band width keeps decreasing.
Remarkably, this second family of flat bands is well-described by an effective p x -p y tightbinding model on a honeycomb lattice, depicted schematically in Fig. 2 (a), and conveniently described by the following Hamiltonian: where the first term describes the σ hopping (head to tail) between the p-orbitals and the second term denotes the π hopping (shoulder to shoulder). Furthermore, n ij = (r i − r j )/|r i − r j |, with r i being the position of site i and n ⊥ ij = U n ij with U being the two-dimensional 90 degree rotation matrix U = 0 −1 1 0 . Finally, t σ and t π ( t N σ and t N π ) are the nearest neighbor (next nearest neighbor) hopping amplitudes for the σ-bonding term and π-bonding term, respectively. (1) upon the choice of hopping parameters t π = 0.25t σ , t N σ = 0.07t σ and t N π = −0.04t σ . The density of states exhibits a characteristic four van Hove singularities structure, with two originating from the Dirac bands and two stemming from the additional two ultra-flat bands. The small ratio between the nearest neighbor hopping amplitudes t π /t σ determines the residual small dispersion in the ultra-flat bands we report. This ratio is controllable by the twist angle, which is summarized in the inset of Fig. 1 (c). All these parameters are related to the interlayer Moiré potential and thus expected to be also affected and controllable by the uniaxial pressure perpendicular to the layers as demonstrated for twisted bilayer graphene (4). We have thus unambiguously established twisted MoS 2 to be a candidate system to realize a p x -p y model on the honeycomb lattice with strongly asymmetric hoppings t σ and t π , giving rise to a new set of ultra-flat bands.

Correlations and magnetic properties
We now study the role of electronic interactions. As the highly-anisotropic p x -p y orbital structure constitutes the essential novelty of twisted bilayer MoS 2 , we focus on quarter filling where orbital fluctuations can be expected to be crucial. This filling fraction is straightforwardly accessible in experiment via back gating, and we defer a discussion of the half-filled case to the Supplementary Material. To proceed, we assume purely local electronic interactions, which can be generically parameterized in terms of the Hubbard-Kanamori Hamiltonian: for two orbitals with rotational symmetry. Furthermore, our DFT calculations suggest t π ≈ 0.25t σ and only weak next-nearest neighbor hopping at small twist angles; we therefore neglect  σ /U , assuming ferro-(blue) or antiferromagnetic (red) order for the spin degrees of freedom. We take the ab initio parameters, t π = 0.25t σ and use an iterative energy minimization. The lower panel determines the phase boundaries for the orbital degrees of freedom given the energetically more favorable spin order shown in the top panel. At J/U = 0.1 we find the spin order to change from AFM to FM, with AFO nematic order for the orbital degrees of freedom remaining stable in agreement with (47). (b) Configurations of orbital vectors found at the end of iterative minimization. Note that we display the projection of τ to the plane in R 3 (indicated by the axis shown in the bottom left), such that nematic states with finite contributions only in xz direction ((1) & (2)) can be distinguished from magnetic states (3) which point perpendicular, i.e along the y-axis.
next-nearest-neighbor hopping in the analysis below. Fig. 3 (c) depicts the local density of states as a function of Hubbard U and Hund's exchange J interactions, calculated via an exact diagonalization study of Eqs. (1) and (2) for a cluster depicted schematically in (a). Clear evidence of a charge gap beyond U/t σ ∼ 4 at small J signifies the onset of a correlated insulator which could be directly observed via transport and scanning tunnelling microscopy. The behavior of the gap is depicted in Fig. 3(b) as a function of U, J and signifies that charge fluctuations are strongly suppressed for large U .
In this regime, a natural follow-up questions concerns possible orderings of the orbital and magnetic degrees of freedom. The corresponding strong-coupling Kugel-Khomskii Hamiltonian (50-52) for the p x -p y model at quarter filling is given in Ref. (45)(46)(47) and reads: Here, ξ 1 ij = 3/4 + S i S j denotes the projector onto triplet states, whereas ξ 0 ij = 1/4 − S i S j selects the singlet spin states instead. Note that the orbital operators, for example Q ij , are bond dependent, giving rise to a strong spatial anisotropy of the resulting spin-orbit model. Details about their definitions can be found in the Methods section.
To study its ground state phase diagram using the ab initio parameters found in the previous section, we employ a mean-field analysis of competing orbital orderings with ferromagnetic and antiferromagnetic spin order. To this end, we note that on the bipartite honeycomb lattice the SU(2) invariant spin sector would, on its own, order either ferro-or antiferromagnetically, depending on the sign of the exchange couplings. As an Ansatz we therefore assume that one of the respective states is stabilized and decouple the spin from the orbital degrees of freedom by replacing S i S j with its expectation value S i S j = ±1/4 such that ξ 1 ij = 1, ξ 0 ij = 0 for ferromagnetic spin order and ξ 1 ij = ξ 0 ij = 1/2 for Neél order.
After such a mean-field decoupling corresponding to the ground state in the spin sector, we analyze the ground states of the resulting Hamiltonian for the orbital degrees of freedom, which we approximate as classical vectors. We use an iterative energy minimization combined with simulated annealing techniques (see Methods) to converge the mean field equations and find the phase diagram summarized in Fig. 4. For example, for our ab initio band structure parameters, a noncollinear spin dimer phase is predicted in a certain range of interaction couplings and even a quantum spin-orbital liquid is found in its proximity (47).

Discussion
We have established that twisted bilayer MoS 2 is a promising platform to realize the orbital anisotropic p x -p y Hubbard model by employing large scale ab initio calculations. We find that

Details on strong-coupling Hamiltonian
Here we give the definitions of the operators used in (3) and discussed in (46). The operators Q ij andQ ij describe processes where orbital occupations of sites i and j are reversed, that is they are defined as with where The orbital projection operators can then be expressed as where e.g. P xx ij selects states where the superposition (p x e x + p y e y )n ij is occupied on nearestneighbor sites i and j connected by the bond n ij .

Details on minimization procedure for classical Hamiltonian
Metropolis Monte Carlo simulations are a prime tool for the investigation of classical spin models, since they allow for off-diagonal, spatially anisotropic spin couplings to be included, even when one-spin terms, such as magnetic fields, are involved. Here we employ a special variant of the algorithm to the mean-field version of (3), keeping in mind that the 'spins' used AFM (spin) Figure 5: Magnetic phase diagram of the strong-coupling Hamiltonian (3) in the isotropic limit t σ = t π . Our results from iterative minimization are in agreement with (46), stabilizing FO nematic order (1) for J/U < 0 (see main text) and FO magnetic order (3) (see main text) for J/U > 1/3. In the intermediate range of parameters 0 < J/U < 1/3 the ferromagnetic spin sector is selected, such that, due to vanishing ξ 0 ij , rotation invariance is restored for the orbital vectors, giving rise to AFM order (2) but without any preferred axis in euclidean space.
in the simulation are approximations to orbital operators τ .
First, a lattice site i is randomly chosen, and its respective gradient field h i = ∇ i H is computed for the current spin configuration {τ i }. Second, a random orientation τ i for the vector at site i is proposed and the weight is computed for an effective inverse temperature β. Performing several Metropolis updates with increasing values of β we are able to efficiently lower the energy of a random initial configuration, minimizing the odds to converge to a local minimum by only allowing optimal updates (i.e. τ i = −h i ) right from the start. After N a sweeps over the full lattice, the so-obtained configuration is ameliorated by N o optimization sweeps, where the randomly selected spin is rotated anti-parallel to the local gradient field such that the energy is deterministically lowered in every step and we converge as close to the global energy minimum as possible. Hence, this algorithm is reminiscent of Monte Carlo simulations with simulated annealing, but at zero temperature where thermal fluctuations are frozen out.
To benchmark our implementation we have carried out the minimization procedure in the isotropic limit t σ = t π for N a = N o = 10 5 , where the optimization sweeps are terminated when the energy change after one sweep, , becomes small (usually ≤ 10 −10 ). Mapping out the phase diagram for both the FM, S i S j = 1/4, as well as the AFM, S i S j = −1/4, spin sector on a lattice with N = 1250 spins subject to periodic boundary conditions we find the result in Fig. 5, which is consistent with the one presented in (46). For J < 0 the AFM spin sector has lower energy, with the orbitals forming a ferro-orbital (FO) nematic state where τ x/z i = 0 and τ y i = 0. For J > 0 one finds the FM spin sector (for which the orbital degrees of freedom restore their rotation invariance) to dominate as long as J < 1/3, where the AFM sector takes over again and establishes a FO magnetic state, i.e. τ x/z i = 0 and τ y i = 0.