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. The origin of these dispersionless bands, is similar to that of the flat bands in the prototypical Lieb or Kagome lattices and co-exists with the general band flattening at small twist angle due to the moiré 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.

T wo-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 another van der Waals materials stacked atop each other at a twist 7-17 displaying features of correlated physics that afford a high level of control. In particular, bi-, tri-, and quadruple-layer graphene 18 as well as twisted few-layer transition metal dichalcogenides (TMDs) 19,20 are currently under intense experimental scrutiny 13,[21][22][23][24][25][26][27][28][29] . 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. Therefore, twisted systems enable the realization of new correlated condensed matter models, establishing a solid-state quantum simulator platform 30 .
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 well-known alternate pathway to flat bands can occur in certain lattices such as the Lieb and the Kagome lattices. Here, purely geometric considerations lead to the formation of perfectly localized electronic states that have weight only on single plaquettes or hexagons, respectively, and that are eigenstates of the kinetic Hamiltonian due to destructive interference between lattice hopping matrix elements 31 . To put it differently, linear combinations of the macroscopically degenerate extended Bloch states in these systems allows to form localized Wannier-like eigenstates (living on single plaquettes or hexagons in the examples above) with no dispersion (for a review on the subject see, e.g, 32 ). 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 32 .
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 33,34 . We confirm that families of flat bands emerge when two sheets of MoS 2 in the 2H structure are stacked at a twist 12,35 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 bandgap 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 36,37 . 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 [38][39][40][41][42] . The strongly asymmetric p x -p y honeycomb model itself 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 in this model originate from destructive interference, in analogy to flat bands in the Lieb and the Kagome lattices 31 discussed above, and will be referred to as ultra-flat bands in the following discussion.
On top of that, the total bandwidth of the strongly asymmetric p x -p y honeycomb effective model realized here (all four bands) can be further flattened by decreasing the twist angle. In addition, these ultra-flat bands can be topologically nontrivial in the presence of spin-orbital coupling (SOC) 43 . Although all the flat bands discussed here originate from the Γ-point states of MoS 2 and are not affected by intrinsic SOC (see Supplementary Fig. 3), we expect that substrate engineering 44 can be used to introduce SOC coupling into these bands and invoke topologically nontrivial behavior of the ultra-flat band states. Previously, the p x -p y model was studied in the context of cold gases where exotic correlated phases were predicted 36,45,46 , as well as in semiconductor microcavities 47 and certain 2D systems such as organometallic frameworks 48,49 and Bismuth deposited on SiC 50 with a focus on their nontrivial topology properties. Our findings elevate tbMoS 2 to an interesting platform where effects of ultraflat 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 whichdepending on the specific parameters and symmetries of the modelhosts magnetic, orbital as well as valence-bond orderings, or even more exotic quantum spin-orbital liquid phases [51][52][53] . With this, our work adds an interesting type of lattice modelthe highly asymmetric p x -p y Hubbard modelto 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 bandwidth of the model (by the twist angle).

Results
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 firstprinciples 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 effects of correlations are expected. The results of such a characterization are summarized in Fig. 1. Note that atomic relaxation has been shown to affect the electronic properties of twisted 2D materials 12,35,54 . While for twisted bilayer graphene this effect is only significant at twist angles smaller than 1 degree 54 , it noticeably alters the low-energy band dispersions and charges density localization for twisted transition metal dichalcogenides bilayer (such as MoS 2 ) even with relatively large twist angles above 1 degree 12,35 . Therefore, we relax all the systems in our DFT calculations. Panel (a) shows the relaxed atomic structure of two sheets of MoS 2 in real space, twisted with respect to each other. A moiré interference pattern forms at a small twist angle yielding a large unit cell, within which we identify different local patterns of stacking of the two sheets of MoS 2 , indicated via areas framed by cyan, magenta or purple dashed lines. The local stacking arrangements of the respective areas are given in the right sub-panels of the panel (a). Note that the B Mo/S and the B S/Mo regions are equivalent as they are related by C2 symmetry. These equivalent B Mo/S /B Mo/S regions form a hexagonal network.
In panel (b) we show the ab initio band structure of the twisted material after relaxation, where we find two families of bands that 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 ∘ -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 bandwidths of the order of 10 meV as the angle approaches Θ ≈ 2°. Similar features are also shown in the work of Naik et al. 35 . The bandwidth and the shape of the flat bands (in particular for the second set) in our calculations are slightly quantitatively different from the previous work probably because we relax the structure directly with DFT while the authors of ref. 35 use a force-field approach. 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 (see Supplementary Fig. 2 for a DFT characterization of the orbital contribution to the different bands). This is different to the case of twisted WSe 2 , where the top valence flat bands originate from the states around the K point in the Brillouin zone of the primitive unit cell (dominated by W d x 2 Ày 2 and d xy orbitals), which experience different interlayer moiré potentials compared with those of the Γ-point flat bands discussed here leading to an effective triangular lattice Hubbard model 13 . Since also 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 55,56 , 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 bandwidth). 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 bandwidth keeps decreasing. To access this second set of bands we need to empty the bands marked in blue first. The effects of this doping are of minor quantitative nature (see Supplementary Fig. 5).
Remarkably, this second family of flat bands is well-described by an effective p x -p y tight-binding model on a honeycomb lattice, depicted schematically in Fig. 2a, and conveniently described by the following Hamiltonian: where c i;s ¼ ðc i;x;s ; c i;y;s Þ T with c i,x(y),s annihilating an electron with p x(y) -orbital at site i and with spin s = ↑, ↓. i; j ( i; j ) denotes (next) nearest neighbors. For each sum in Eq. (1), 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 k ij ¼ ðr i À r j Þ=jr i À r j j, with r i being the position of site i and n ? ij ¼ Un k ij with U being the two-dimensional 90 degree . 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. Figure 2b, c depict the corresponding dispersions, density of states, and wave functions in comparison to model predictions, illustrating that the four moiré bands at low energies are well captured by Eq. (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. 1c. All these parameters are related to the interlayer moiré potential and are thus expected to be also affected and controllable by the uniaxial pressure perpendicular to the layers as demonstrated for twisted bilayer graphene 4 .
The flat band wavefunctions consist of atomic wavefunctions from the p z orbital on S atoms and the d z 2 orbital on Mo atoms. Modulated by the moiré potential, the weighting of the atomic wavefunctions and their modulus square (i.e., charge density) vary at different atomic sites across the whole supercell, showing distinct patterns for different flat band states at the K point in the supercell Brillouin zone as shown in Panel (c) of Fig. 2. These patterns of the charge density as well as the real and the imaginary part of the total wavefunctions obtained from DFT show features consistent with those of the p x -p y Hamiltonian of Eq. (1). Note, that we call this the p x -p y Hamiltonian to connect to established literature on the subject; whereas the actual moiré wave functions are composed of p z and d z orbitals, they transform like p x , p y orbitals according to the irreps of the reduced symmetry group of the moiré supercell. Interestingly, the charge density distribution of the top ultra-flat band state displays a Kagome lattice structure. 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 (one electron per sublattice in the Moié unit cell) where orbital fluctuations can be expected to be crucial. This filling fraction is straightforwardly accessible in the experiment via back gating, and we defer a discussion of the half-filled case to Supplementary Note 1. 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. More realistic modelling should include long-range interactions. However, for our choice of commensurate quarter filling, any longer-ranged component of the Coulomb interaction at strong-coupling will serve merely to renormalize the effective spin-orbital interactions of the resulting Kugel-Khomskii model and we therefore concentrate on purely local interactions for simplicity. Furthermore, our DFT calculations suggest t π ≈ 0.25t σ and only weak next-nearest neighbor hopping at small twist angles; we therefore neglect next-nearest neighbor hopping in the analysis below (see Supplementary Fig. 4 for a comparison of the band structures with and without next-nearest neighbor hopping). An ab initio based characterization of the values of U and J requires numerically expansive Wannierzation of the wave functions and is unfortunately beyond the scope of this work. However, by substrate engineering 22 it is likely that a whole range of values can be accessed and therefore it is useful to vary these parameters to explore all possible phases accessible in experiments to make concrete predictions. Vice versa given a future experimental observation our results can be used to estimate the strength of correlations. Figure 3c 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. 3b as a function of U, J and signifies that charge fluctuations are strongly suppressed for large U. Establishing the existence of a charge gap motivates to set up a strong-coupling Hamiltonian routinely employed for the types of systems under scrutiny here.
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 57-59 for the p x -p y model at quarter filling is given in refs. 51-53 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. To be more precise following ref. 52 , the operators Q ij and Q ij describe processes where orbital occupations of sites i and j are reversed, that is they are defined as The orbital projection operators can then be expressed as P xx , where e.g. P xx ij selects states where the superposition ðp x e x þ p y e y Þn k ij is occupied on nearest neighbor sites i and j connected by the bond n k ij . 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 for orbital orderings with ferromagnetic and antiferromagnetic spin order. Note, that the simplifying assumption of vanishing temperaturea standard one in condensed matter researchstill allows to draw conclusions for the low-temperature physics accessible in experiments as fingerprints of the phases we discuss extend into this regime as well. 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. Panel (a) shows the energy of ferromagnetic and antiferromagnetic spin configurations from which the magnetic phase diagram can be read off. This is given in the upper part of the plot and we find antiferromagnetic ordering with an intermittent ferromagnetic phase at intermediate ratios of 0.1 < J/U < 1/3. In the lower part of the plot, we show the corresponding subsidiary orbital order. From our simulations, we identify three different configurations of orbital vectors τ, which can be classified according to their projection on a single definite plane in space, shown in the lower left of the plots: (1) ferro-orbital (FO) nematic order 5,6,60-62 , where the vectors on all lattice sites align in parallel to the xzplane. Quantum mechanically, finite values of hτ x=z i i indicate an imbalance of the occupation of p x and p y orbitals, breaking rotation symmetry and thereby motivating the notion of a nematic state. (2) AFO nematic order; each vector is aligned antiparallel with its nearest neighbors corresponding to hτ x=z i i ≠ 0 on each sublattice, but without finite projections τ y i on individual sites. (3) FO magnetic order; all vectors order along the y-axis, such that hτ y i i ≠ 0, which, in the quantum mechanical system, would indicate time-reversal symmetry breaking. The inclusion of quantum fluctuations can change this picture and more exotic ground states may emerge. 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 53 . Since these exotic phases primarily occur for weak Hund's coupling and strong orbital anisotropies, the assumptions made for our calculations can therefore be justified for sizable J H and modest distances to the isotropic t σ = t π point.

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 families of flat bands emerge where the first family of flat bands shows s-orbital character and the second family is an intriguing realization of a strongly asymmetric p x -p y Hubbard model both on a honeycomb lattice, adding a lattice with nontrivial almost perfectly-flat bands due to destructive interference to the growing list of systems that can be engineered in twisted heterostructures. The symmetry of these flat bands is inherited from the hexagonal lattice formed by the equivalent B Mo/S and B S/Mo regions. At an even smaller angle, the sequence in the family of flat bands found with respect to their orbital character continues. Our analysis shows that the low-energy DFT band structures in this system can σ =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 ref. 53 . b Configurations of orbital vectors are found at the end of iterative minimization. Note that we display the projection of τ to the plane in be well captured by a free electron gas model modulated by a simple harmonic potential that has hexagonal (D 6 ) symmetry, which is consistent with a recent study 63 . This simple model further shows that the next family would exhibit a d-orbital character on the honeycomb lattice. Such a lattice would effectively realize a multi-orbital generalization of a Kagome latticea prototypical model for quantum spin liquids. However, at such small angles strong relaxation is likely to become dominant, prohibiting access to this regime and potentially spoiling its experimental realization. Currently, the ab initio characterization of such small angles is numerically too exhaustive and this work sparks a direct need for novel computational methods to tackle this question.
Furthermore, our combined exact diagonalization and strongcoupling expansion approaches classify the magnetic and orbital phase diagrams, however, the inclusion of quantum fluctuations stipulates an intriguing avenue of future theoretical research.
Indeed, previous theoretical works provide some evidence for a quantum spin liquid in the SU(4)-symmetric Kugel-Khomskii model on the honeycomb lattice 64 , the square lattice as a related system without frustration 65 and studied the role of perturbations that break SU(4) symmetry and isotropy 53 In twisted MoS 2 , this regime would in fact map to larger twist angles, where the anisotropy of the p x -p y model is less pronounced, as well as to a regime of vanishing Hund's coupling, placing such a putative quantum spin liquid at the transition between FO nematic and AFO nematic phases.
In addition, by proximity or variations in the chemical composition of the twisted bilayer, it might be possible to induce spinorbit coupling splitting of the ultra-flat bands at the top and bottom of the asymmetric p x -p y dispersion. Such a bandgap opening would induce interesting topological properties 66 in a highly tunable materials setting.

Methods
Details on ab initio calculations. We calculate the electronic properties of twisted bilayer MoS 2 with ab initio methods based on density functional theory (DFT) as implemented in the Vienna ab initio Simulation Package (VASP) 67 . We employ plane-wave basis sets with an energy cutoff of 550 eV and pseudopotentials as constructed with the projector augmented wave (PAW) method 68 . The exchangecorrelation functionals are treated at the generalized gradient approximations (GGA) level 69 . The supercell lattice constants are chosen such that they correspond to 3.161 Å for the 1 × 1 primitive cell of MoS 2 . Vacuum spacing larger than 15 Å is introduced to avoid artificial interaction between the periodic images along the z-direction. Because of the large supercells, a 1 × 1 × 1 k-grid is employed for the ground state and the relaxation calculations. For all the calculations, all the atoms are relaxed until the force on each atom is less than 0.01 eV/Å. Van der Waals corrections are considered with the method of Tkatchenko and Scheffler 70 . We extract the real and the imaginary parts of the DFT wavefunctions with the VASPKIT code 71 .
Details on exact diagonalization. Exact diagonalization calculations were performed for the electronic tight-binding model in Eq. (1) with Hubbard-Kanamori interactions defined in Eq. (2). All calculations were performed for a two-orbital eight-site cluster with periodic boundary conditions at quarter filling, corresponding to eight spin-1/2 particles in sixteen orbitals. Rotationally symmetric Kanamori interactions are adopted, with U 0 ¼ U À 2J. As the magnitudes of the Hubbard U and Hund's exchange J interactions cannot be reliably predicted for a Moié supercell from first principles, all presented results are shown as a function of U, J. Calculations of the single-particle Green's functions and local density of states are performed starting from the ground state in the total momentum K tot = 0 and total spin S z = 0 sectors, using the Lanczos method and continued-fraction representation, and a spectral broadening (imaginary part of the self-energy) of η = 0.1 is imposed.
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 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 τ 0 i for the vector at site i is proposed and the weight g ¼ min e Àβðτ 0 i Àτ i Þh i ; 1 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 ref. 52 . For J < 0 the AFM spin sector has lower energy, with the orbitals forming a ferro-orbital (FO) nematic state where hτ x=z i i ≠ 0 and hτ y i 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. hτ

Data availability
The raw data sets used for the presented analysis within the current study are available from the corresponding authors on reasonable request.

Code availability
The tailored developed codes used in this work can be provided from the corresponding author on reasonable request. Ab initio calculations are done with the code VASP (version 5.4.4).
Received: 10 December 2020; Accepted: 1 September 2021; Fig. 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 ref. 52 , 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.