The Orbitally Selective Mott Phase in Electron Doped Twisted TMDs: A Possible Realization of the Kondo Lattice Model

Moir\'e super-potentials in two-dimensional materials allow unprecedented control of the ratio between kinetic and interaction energy. By this, they pave the way to study a wide variety of strongly correlated physics under a new light. In particular, the transition metal dichalcogenides (TMDs) are promising candidate"quantum simulators"of the Hubbard model on a triangular lattice. Indeed, Mott and generalized Wigner crystals have been observed in such devices. Here we theoretically propose to extend this model into the multi-orbital regime by focusing on electron doped systems at filling higher than 2. As opposed to hole bands, the electronic bands in TMD materials include two, nearly degenerate species, which can be viewed as two orbitals with different effective mass and binding energy. Using realistic band-structure parameters and a slave-rotor mean-field theory, we find that an orbitally selective Mott (OSM) phase can be stabilized over a wide range of fillings, where one band is locked in a commensurate Mott state, while the other remains itinerant with variable density. This scenario thus, realizes the basic ingredients in the Kondo lattice model: A periodic lattice of localized magnetic moments interacting with metallic states. We also discuss possible experimental signatures of the OSM state.

In contrast, semiconducting transition metal dichalcogenides (TMDs) subject to moiré potentials are expected to have a simpler microscopic picture. The low energy physics is captured by a Hubbard model on a triangular lattice 19,20 . Mott insulators and generalized Wigner crystals have been experimentally observed [21][22][23][24][25] , as well as possible indications of superconductivity 21 . The relative simplicity of their microscopic starting point makes the TMD moiré devices prime candidates for condensed matter "quantum simulators" of the Hubbard model.
A canonical model that is both of great fundamental interest to quantum condensed matter, and has not yet been realized in moiré devices, is the Kondo lattice model 26 . Its main ingredients are a lattice of localized moments coupled to a Fermi liquid of itinerant electrons. The main coupling between these two degrees of freedom is spin-exchange. The case where the strongest exchange mechanism is antiferromagnetic is understood to be the minimal model that captures the low-energy physics of many rare-earth compounds, known as heavyfermion materials 27,28 . When the dominant exchange is the Hund's coupling between the local and itinerant orbitals, the coupling is ferromagnetic. Such a scenario was discussed in the context of the orbitally selective Mott (OSM) phase 29,30 . Materials that host coexisting itinerant and localized states, exhibit a plethora of exotic phases such as heavyfermi liquids, metallic magnets, high-T c superconductors and non-Fermi liquids [29][30][31][32][33][34][35][36][37][38][39][40][41] . However, what makes them arXiv:2103.06313v3 [cond-mat.str-el] 26 Oct 2021 species a species b E [eV] ' FIG. 2. Left: Typical dispersion of the two lowest Bloch bands of the conduction band in a bare single-layer TMD obtained from the tight-binding model 45 . Near the K and K points the bands are approximately parabolic and assume a small splitting due to spin orbit coupling in the second order (see inset). Right: Upon lightly doping the system two Fermi pockets of opposite spin orientation form around each high symmetry point, corresponding to the two Bloch bands.  31 36 especially interesting is the existence of quantum phase transitions, where the lattice of local moments melts into a metallic state 25,31,34,37,40,[42][43][44] . Such a transition is not captured by the Ginzburg-Landau paradigm because it must include a whole Fermi surface that emerges at the quantum critical point 35,38,41 . What controls the different ground states, and the nature of the quantum critical point separating them, is still debated. However, the comparison between theory and experiment becomes challenging due to the complex structure of the materials which realize this physics. For this reason, a controlled experimental realization of such a minimal model is highly desirable. In this paper we explore the conditions under which, TMDs subject to a moiré potential can host a state of coexisting itinerant and localized electrons known as the orbitally selective Mott state (OSM). 30,41,46 We first argue that the mini-band structure of electron doped TMD moiré devices can potentially host multiple flat bands, which can be simultaneously at a state of partial filling. This is mainly because of the relatively small spinorbit splitting of the bare conduction bands around the K and K points 45 [see Fig. 1.(a)]. We consider the situation where such minibands are induced in one layer by another "inactive" layer in a heterogeneous structure 19 . For example, two prototypical bilayers we consider are WS 2 /WSe 2 , where the effects of spin-orbit splitting on the conduction bands are small but noticeable, and MoS 2 /MoSe 2 , where the splitting is negligible. In FIG. 3. The two lowest minibands of the two species near the conduction band minima resulting from a moiré potential of angle θM = 3 • and depth V0 = 15 meV. The dashed lines indicate the average band energies,¯ τ (see text). In this paper we focus on the lowest pair of such minibands that overlap in energy space. (a) For WS2/WSe2 this pair includes the lowest miniband of species b and the first excited band of species a.
(b) For MoS2/MoSe2 the SOC is very weak and the miniband structure is nearly degenerate for all minibands. Therefore, for this material we consider the lowest miniband for both Bloch bands.
both cases the sulfur based compounds are where the electronic states reside, and the selenium based layers take the role of the "inactive" layer that induces the moiré potential. Using a slave-rotor mean-field approximation 47-49 within a simplified on-site interaction Hubbard model, we identify the emergence of the OSM state at fillings surrounding n = 4 or n = 2 (depending on the strength of spin-orbit coupling). In this phase one species is in a Mott state, while the other species is partially filling one of its minibands [see

II. Model Hamiltonian
The two lowest Bloch bands above the band gap, which we denote here as species τ = a, b, have nearly degenerate band minima in the vicinity of each valley, K and K 45 (see Fig. 2). Due to spin-orbit coupling they are split and assume different effective masses. As mentioned above, this splitting is significantly smaller compared with the equivalent splitting in the valence band. Nonetheless, the spin projection along z remains a good quantum number up to second order in perturbation theory.
As usual we obtain an additional valley degree of freedom by expanding the momentum around K and K . Note however, that spin-orbit coupling slaves spin to valley within a given Bloch band. We therefore denote the additional degree of freedom by its spin as follows: For the lower Bloch band τ = a, the state σ =↑ and σ =↓ corresponds to a valley K and K , respectively. On the other hand, for the higher Bloch band τ = b, the state σ =↑ and σ =↓ corresponds to a valley K and K, respectively.
The resulting Hamiltonian (up to quadratic order in deviation from the high-symmetry points), is given by Where ∆ τ and m τ are the species dependent band minimum and mass, respectively. The values are listed in Table I. In principle, the Fermi surfaces surrounding the K and K points (in terms of the momentum relative to these points), are non-degenerate except for six high symmetry lines. However, as a result of the parabolic band approximation, used in Eq. (1), the Fermi surfaces are spherically symmetric and are thus doubly degenerate everywhere (the degeneracy corresponds to the spin index σ). This reflects an emergent SU(2) symmetry of each species 19 .
In this paper we consider two prototypical TMD bilayers WS 2 /WSe 2 and MoS 2 /MoSe 2 . In both cases the band alignment properties are such, that the charge carriers reside on the sulfur based side of the bilayer upon electronic doping. Thus, the selenium based layers are inactive and only induce the moiré potential. In the case of MoS 2 /MoSe 2 spin-orbit coupling is very weak and consequently the masses m τ and band minimum points ∆ τ in Eq. (1) are almost identical. In the case of WS 2 /WSe 2 , the effects of spin-orbit coupling are more noticeable, such that m b /m a ≈ 0.75 and ∆ b − ∆ a = 30 meV.
We now turn to consider the influence of a moiré potential on the band structure close to the bottom of the conduction bands τ = a, b. We follow Refs. 19 and 20. The induced potential is given by where G G G j =R(j π 3 )(4π/ √ 3a Mx ), j = 0, . . . , 5 are the six shortest reciprocal lattice vectors of the moiré super lattice. a M = a/θ M is the moiré lattice constant and θ M = √ δ 2 + θ 2 is the effective twist angle. Here δ is the lattice mismatch and θ accounts for any additional twist. We take the strength of the potential to be V 0 = 15 meV for both bilayers.
The parabolic Hamiltonian Eq. (1) together with the moiré potential Eq. (2) are diagonalized using a nearlyfree electron approximation truncated at the level of 19 bands (3rd nearest neighbour in reciprocal space).
In Fig. 3 we plot the two lowest minibands of each species using realistic parameters for the bilayers. In panel (a) we show that for the strongly spin-orbit coupled bilayer, WS 2 /WSe 2 , the lowest miniband of species b overlaps with the first remote miniband of species a. On the other hand, in panel (b) we show that for the weakly spin-orbit coupled bilayer, MoS 2 /MoSe 2 , the minibands of the two species are almost identical. In this case the two lowest minibands (and the two first excited bands) overlap in energy.
We will be interested in the physics arising from partially filling two different minibands simultaneously. Therefore, from here on we will focus exclusively on the lowest pair of mninbands that overlap in energy space corresponding to the two Bloch bands τ = a, b. The miniband Hamiltonian then assumes the form However, we still define the density in units of total filling starting from the bottom of the conduction band. Consequently, for WS 2 /WSe 2 [ Fig. 2(a)] the relevant range of filling is n ∈ [2,6], where the lowest miniband of species a is already completely filled and contributes a background charge of 2. This situation is also depicted in the center of panel (b) in Fig. 1. On the other hand for MoS 2 /MoSe 2 the two lowest minibands of each species overlap and therefore we focus on the range of filling n ∈ [0, 4]. The second ingredient in our model is the interaction. We consider an on-site repulsion of the form where δn τ = σ ψ † τ σ ψ τ σ − 1 is the density operator in particle-hole symmetric form. 50 η is a phenomenological parameter, which accounts for the possible difference in the Wannier-orbital spread of the species. When the lowest miniband of species b overlaps with the first remote miniband of a ( Fig. 3. a), the spread of the Wannier orbital of the latter is expected to be wider than that of the former. Consequently, electrons in miniband a will have a weaker Coulomb repulsion, corresponding to η < 1. On the other hand, when the overlapping minibands are both the first flat band ( Fig. 3. b) the interactions are expected to be roughly equal and η = 1. We consider a constant interaction U of moderate strength, which corresponds to the estimate of Ref. 19 with a large dielectric constant κ ≈ 5 51 . 52

III. Slave Rotor Mean-Field analysis
We now turn to study the ground state of the Hamiltonian Eqs. (3,4). In particular, we are interested to understand whether a Mott state of one of the species can be stabilized over a finite density range, where the other band remains metallic. To this end, we employ the slaverotor mean field theory [47][48][49] . It consists of decomposing the field operators into bosonic rotors multiplied by neutral spinon operators ψ iτ σ = e −iθiτ f iτ σ . The respective density operator of each species, which are conjugates of the phases above, are then represented by angular momentum operatorsL iτ = −i∂/∂θ iτ , subject to the local where the sum over spin is implicit.
The above decomposition allows for a mean-field treatment of the Mott transition 47 . The corresponding "order parameter" is the quasiparticle weight Z τ = | e iθτ | 2 . When the rotor fields are pinned Z τ = 0, resulting in a finite overlap between the quasiparticle and bare-electron operators. Moreover, the uncertainty principle implies the conjugate charge operatorL τ experiences large fluctuations. Thus, we can identify this phase with a Fermi liquid. On the other hand, when the charge operatorŝ L τ are pinned, which corresponds to small charge fluctuations, the conjugate phases are strongly fluctuating and Z τ = 0. This phase is thus associated with the Mott state.
Before applying the slave-rotor decomposition however, it is essential to decompose the miniband dispersion, ξ kτ , Eq. (3) into two terms where¯ τ is the average energy of the miniband (¯ τ = k∈M BZ ξ k,τ ), and the remainder, kτ , is the kinetic part of the dispersion, which averages to zero.¯ τ can be interpreted as the effective binding energies of electrons to the respective minibands (dashed lines in Fig. 3). The importance of this decomposition, is to separate these local energy shifts from the dispersion because they should not be renormalized by the quasiparticle weights Z τ . Indeed, in the slave-rotor theory the quasiparticle weight only renormalizes the band width but does not shift the average energy of the band 47 .
Performing the slave-rotor decomposition to both species we obtain the Hamiltonian where t ij τ are the set of tight-binding parameters that reproduce the dispersive part k in Eq. (2) when transformed to reciprocal space.
To asses the ground state of the Hamiltonian Eq. (6) we employ the variational method, as opposed to Ref. 47, where the self-consistent mean-field approach was used. Namely, we minimize the expectation value of Eq. (6) with respect to the variational wavefunction denoted by This variational state is a product of the ground-states of the rotor Hamiltonians and two Slater-determinant states ("Fermi sea" states) of the spinons, where the density is controlled by the chemical potentials µ a and µ b . We must determine six variational parameters with three constraints L τ = f † iτ f iτ −1 and τ f † iτ f iτ = n. The parameter K τ controls whether species τ is metallic or localized. When the minimal energy solution is obtained with K τ = 0 the rotors are pinned and we get a finite quasiparticle weight Z τ = 0 corresponding to the metallic state. On the other hand, for K τ = 0 the rotors are in eigenstates of the angular momentum operator where the average of e iθτ vanishes, corresponding to the Mott state (Z τ = 0). Additionally, there is a freedom to redistribute charge between the two bands, which is controlled by the difference in the chemical potentials, µ a − µ b . Finally, the constraints are fulfilled using the three Lagrange multipliers h τ and the sum µ a + µ b . Note that in the case of WS 2 /WSe 2 , where the overlapping bands include one the first excited band of species a, which is more disprersive compared to the lowest band of species b, we apply the slave-rotor decomposition only to band b. For more details on the slave-rotor analysis we perform and the minimization procedure see Appendix C.
In Fig. 4 we plot the phase diagrams resulting from the variational minimization. Panels (a)-(d) correspond to the WS 2 /WSe 2 bilayer using U = 60 meV and η = 1/2. Panel (a) and (b) are maps of the filling of species b and a, respectively, in the space of total filling n and twist angle θ M . In this case we recall that there is another completely filled miniband below the relevant pair of overlapping minibands (see Fig. 3.a), therefore the total filling is given by n = 2 + n a + n b . We turn our focus to the region inside the white dashed line, where the filling of species b is locked to unity, while the filling of species a varies continuously. 53 In the same region, we find that the quasiparticle weight Thus, this region is identified as the OSM phase, where a lattice of localized magnetic moments coexists with itinerant electrons. Panel (c) shows the filling of each band for a specific twist angle θ M = 3 • showing that the density of the itinerant band can be tuned over a large range inside the OSM phase.
In panels (e)-(h) we plot the results for MoS 2 /MoSe 2 using U = 40 meV and η = 1. Panels (e) and (f) are maps of the filling of species a and b, respectively. Note that in this case the total density is simply n = n a + n b Panels (g) and (h) are the quasiparticle weights of bands a and b, respectively. Here we identify two OSM phases, one where band a is locked in a Mott state (OSM a ) and one where band b is locked in the Mott state (OSM b ). Additionally, at n = 2 there is a Mott state of both species, which is expected to have an approximate SU(4) symmetry [54][55][56][57] .
For both materials, the OSM state assumes a large portion of the phase space (in Appendix C we show that this scenario is relevant to other TMD materials). For small twist angles, the density of itinerant electrons can be tuned between completely empty and completely filled states, which potentially allows to tune the strength and sign of the RKKY interaction between local moments, which is mediated by the itinerant miniband 58 . At larger twist angles (θ M 4 • ) the lattice of localized electrons melts into a Fermi liquid. Such a transition is characterized by the emergence of Fermi surface, which is not captured by the Ginzburg-Landau paradigm and is therefore of special interest 35,38,41 . 59 We have also tested the stability of the OSM phase to variations in the parameters ∆¯ =¯ a −¯ b and η numerically. In Appendix D we show that the range of filling where the OSM phase occurs is large for a wide range of ∆¯ and η. We also roughly estimate this range analytically and find that it is expected to be large in the parameter regime |∆¯ | < U (1 − η) 2 /2 − T it . Here T it < 0 is the the kinetic energy gain of filling the Fermi sea of the itinerant band minus the interaction energy associated with onsite fluctuations of charge. This analysis shows that the existence of a wide OSM phase is robust to the parameters of our model.

IV.
Experimental consequences of the OSM state We turn to discus experimental consequences of the OSM phase. We first discuss the enlarged entropy associated with the formation of local moments. When the moments are free they contribute one k B per lattice site. If a magnetic ordering is present, the local-moment contribution will be significant above the ordering temperature. Additionally, a distinct feature of this contribution will be a strong dependence on magnetic field. Indeed, the authors of Refs. 60 and 61 have recently measured such an enlarged entropy in TBG, where they attributed it to local moments coexisting with metallic states. Similarly, in the regime where both phases are metallic, but close to the OSM regime we may expect a Pomeranchuk effect upon heating 60,61 .
To estimate the change in entropy across the OSM transition we assume the local moments contribute their maximal entropy, while the metallic states contribute s where N kτ are the momentum space Fermi-Dirac distribution functions, which include the effects of the quasiparticle weight Z τ . In Fig. 5   per moiré lattice site as a function of density for the WS 2 /WSe 2 bilayer at θ M = 2.5 • for three different temperatures. The distinct signature is a large jump at the boundaries of the OSM state, where we also observe an enhanced specific heat manifested in the strong dependence of S with T .
Another suitable probe for the OSM state is magnetotransport 41 , especially given that we predict this state at a relatively high density range n ∼ 2 − 4, where the effects of disorder are less prominent as compared with the filling range of the lowest miniband. Inside the OSM phase the Hall number, which is seen both in the slope of the classical Hall resistivity and in the period of quantum oscillations, will correspond to a "small" Fermi surface (of volume n b = n − 1). At the phase transition point [green hue in Fig. 1.(c)] the local moments melt into a metallic state, manifested in a Lifshitz transition, where we can distinguish two scenarios. When the dominant exchange interaction between the two species is antiferromagnetic we expect a heavy-fermi liquid state to emerge between the fully metallic phase and the magnetic metal. In this case the Hall number changes from the "small" volume n − 1 to the "large" volume n. The second scenario, is where the exchange interaction is dominated by ferromagnetic exchange (e.g. due to the orbital Hund's coupling). In this case, theory does not predict the emergence of a hybridization gap between the local and itinerant electrons. Instead, a new Fermi surface emerges at the transition point. Thus, we expect the appearance of beating in quantum oscillations and non-linearity of the classical Hall effect (see for example Ref. 62). Thus, magneto-transport measurements across the melting transition can also distinguish the nature of magnetic exchange mechanism.

V. Summary
We proposed that electron doped TMDs subject to a moiré potential are prime candidates to realize the Kondo lattice model. The essential ingredient is the multiplicity of electron band minima close to the K and K points, which allows for two moiré bands of different width to be simultaneously at partial filling. We used a simplified model with constant on-site Coulomb repulsion and a slave rotor mean-field theory to study the possible ground states of the system. We found a large phase space, where an orbitally selective Mott phase forms. Such a state is characterized by a Mott state of one species coexisting with a metallic state of the other. This opens a path to simulate the Kondo lattice model and possible exotic phase transitions in TMD moré devices.
Note added -Upon completion of this paper we came to learn about a related theoretical proposal regarding tri-layers of twisted graphene sheets 63 .

VI. Acknowledgments
We are grateful to Erez Berg, Debanjan Chowdhury, Rafael Fernandes, Efrat Shimshoni, Inti Sodemann and Arun Parameknti for helpful discussions. This research was funded by the Israeli Science Foundation under grant number 994/19. JR acknowledges the support of the Alon fellowship awarded by the Israel higher education council.

A. Continuum dispersion
In this appendix we provide additional information about the computation of the continuum Hamiltonian. We start with the tight-binding approximation for single layer semiconducting TMDs of the trigonal prismatic structure (H) 45 . This model consists of three orbitals d z 2 , d xy and d x 2 −y 2 , taking into account spin-orbit coupling and hopping up to the third nearest neighbours on the triangular lattice. The conduction band consists of two Bloch bands denoted by τ = a, b, which are plotted in Fig. 2 (colored red and blue, respectively) and will be referred to as "species" henceforth. Each such band has two parabolic minima near the K and K corresponding to spin states σ =↑↓ (valley and spin are locked. However, it is important to note that the spin orientations near K and K are opposite in the two Bloch bands). Up to quadratic order in deviations from the high-symmetry points we obtain the Hamiltonian Here ∆ τ and m τ are species dependent band minimum and mass, respectively. k is the lattice momentum relative to the high symmetry points, i.e. relative to K for (a, ↑), (b, ↓) and relative to K for (b, ↑), (a, ↓). A crucial feature, which is unique to the conduction bands, is that the higher order spin-orbit splitting, |∆ a − ∆ b |, is comparable to the expected moiré lattice depth and resulting miniband width (see Table. I).
We now turn to consider the effect of a moiré potential, which we assume is induced by a second layer. At small twist angles θ M π the superlattice constant is given by a M ≈ a/θ M , where θ M ≡ √ δ 2 + θ 2 , δ is the miss-match between the layers taken from Ref. 20 and θ accounts for any additional twist. In this limit we have a M a, which justifies the use of a simple triangular periodic potential constructed out of the six smallest harmonics The potential has three fold rotational symmetry which states that: To obtain the miniband structure of the lowest mini-bands we use a 19 band model without counting degeneracy of spin and species. For simplicity we take the moiré potential strength to be uniform across platforms and given by V 0 = 15 meV 19,20 .
In Fig. 3 we compare the two lowest minibands for the two species (species a colored blue, and species b colored red) for realistic parameters of two candidate materials. As can be seen a feature of these miniband structures is the overlap of bands belonging to different species. Note that the overlapping minibands are not necessarily the same numeral sub band. As shown in panel (a) for WS 2 the overlap is between the first excited band of one species and the lowest miniband of the other (the same is true for WSe 2 and MoSe 2 ). On the other hand for MoS 2 the overlap occurs between the lowest minibands of the two species (panel b). Upon restriction to the two bands of interest (namely, those that are overlapping) we obtain the dispersion in Eq.(3).

B. Interactions
In the paper we assume a contact interaction of the form where δn iτ = ψ † iτ ψ iτ −1. Notice that we have written the interaction in a particle-hole symmetric manner, which can be absorbed into the parameters¯ τ in Eq.(5).
The relative strength of the interaction parameters U τ τ depend on the spread of the Wannier orbitals of the corresponding minibands 19,20 . Namely, when both minibands are the lowest sub-band of their corresponding species [as shown for MoS 2 in Fig.3 (b)], the spread of the two Wannier functions is approximately the same, and we expect U aa U ab U bb . In this case the interaction (B1) is proportional to the square of total density.
On the other hand, when the two overlapping bands belong to different sub-bands [see WS 2 , in Fig.3 (a)] their corresponding Wannier functions will differ in width (namely, the higher, more dispersive band, will have a larger spread). Thus in this case, the interaction parameters may differ significantly. To account for this effect we consider the phenomenological parameter η such that U bb = ηU ab = η 2 U aa . The interaction Eq. (B1) then assumes the form η < 1 describes the scenario where the Wannier function of band b has a smaller spread when compared to a. The value of U itself is twist-angle dependent 19 . For simplicity however, we will take a constant value U = 60 meV for WS 2 , WSe 2 and MoSe 2 , which corresponds to a dielectric environment of ε = 5 51 . For MoS 2 we use U = 40 meV. We note these values are weaker than those used in other studies estimates 20,64 .
The quadratic form of the interaction (B2) was chosen for simplicity. In general, the ratio between the interand intra-species interactions is not controlled by a single parameter η. Therefore, it is important to note that the OSM phase space is expected to be reduced in the case where the interspecies interaction U ab becomes much larger than U aa or U bb . As we show in Appendix D and in the analysis of MoS 2 with =1however, the OSM state is not very sensitive to large interspecies interaction. Another crucial interaction we have neglected is longer range interaction. We expect these interactions to cause additional "Wigner crystal" insulating phases to appear. They will likely cause the phase space of the OSM state to shrink as well. However, these incompressible states may also stabilize over a finite range of doping with the aid of a background incompressible state, i.e. forming an orbitally selective wigner crystal. Finally, we have also neglected spin exchange interactions (e.g. Hund's), which will be discussed in Appendix E.

C. Details of the variational minimization of the slave-rotor mean-field free energy
In this section we describe in more detail the slaverotor mean-field theory [47][48][49] , that we have used in the main text. We first decompose the field operators into bosonic rotors (e −iθiτ ) multiplied by neutral spinon operators (f iτ ) The respective density operator of each species, which are conjugates of the phases above, are then written in terms of angular momentum operatorŝ subject to the local constraint where the sum over spin is implicit. The application of the slave-rotor decomposition to Eq.(3) and Eq.(4) of the main text enables a simple mean-field analysis, which captures the localization-delocalization transition of a half-filled band. This decomposition allows for a mean-field treatment of the Mott transition. The "order parameter" is the quasi-particle weight Z τ = | e iθτ | 2 . When the rotor's phase θ τ assumes a finite expectation value, Z τ = 0 and the spinon quasi-particles have a finite overlap with the original electronic operator. This phase thus, corresponds to a Fermi-liquid state. On the other hand, when θ is delocalized the rotor e iθτ has a vanishing expectation value and the quasi-particle weight disappears (Z τ = 0), corresponding to a Mott phase.
Below we will describe two approaches. First, we will consider the case where only one species is decomposed (the flatter of the two). This case is more applicable to situation, where the density of states of the two bands differ significantly. In the second case, we will consider the same analysis where both bands are decomposed.

Slave rotor decomposition of a single band in a two band system
Applying the aforementioned slave-rotor decomposition to the flatter band (for the purpose of the discussion let it be, τ = b), as is the case for Eq.(3) and Eq.(4) we obtain where t ij τ are the set of tight-binding parameters that reproduce the dispersive part k in Eq.(5) when transformed to reciprocal space.
To estimate the location of possible Mott phases of Eq.(C3) we use a variational approach. This should be cotrasted with Ref. 47 where the self-consistent meanfield technique was used. In the variational approach we minimize the expectation value of Eq. (C3) with respect to a variational wave function denoted by |Ω V = |K b , h b ⊗ |µ b ⊗ |µ a , which is a product of the groundstates of the following variational Hamiltonians where we have shifted the energies such that the center of band b is at zero and ∆¯ =¯ a −¯ b . Eq.(C4) controls the rotor field, where the term proportional to K b acts to pin the phase θ b giving rise to a finite quasi-particle weight Z b . Thus, we can identify the Fermi-liquid (Mott) phases with situations where the minimal energy solution is obtained with K b = 0 (K b = 0). The parameter h b is used to obey the slave-rotor constraint on average. The second and third variational Hamiltonians Eq. (C5) and (C6) generate Fermi sea states of spinons and b-electrons, with density controlled by the parameter µ a and µ b , respectively. Notice that ground state of Eq. (C5) is independent of the band width and therefore Z b is omitted. We then minimize the expectation value of the full Hamiltonian Eq.(6), denoted by with respect to the four parameters K b , h b and µ b and µ a subject to two constraints The difference between the number of constraints and variational parameters implies that two are free. These correspond to the the quasi-particle weight of band b and any distribution of the total density between the bands. These two parameters are dictated by energetics. Notice that in using Eq. (7) we have neglected spatial fluctuations of the field θ b . This restricts our ground state manifold (for example, it can not capture spincorrelations 48 ). However, it allows for a significant simplification: The expectation value of the rotor correlation function becomes a product of local expectation values e i(θiτ −θjτ ) = e iθiτ e −iθjτ = Z τ . Consequently, the expectation value of the kinetic energy terms can be straightforwardly transformed back to momentum space, reproducing the exact continuum dispersion relation Eq.(3).
and N 0 (x) = 1/(e βx + 1). Here β is the inverse temperature, which will be taken to infinity β → ∞, which is used as a numerical parameter to smoothen the discretization.
In panels (a)-(d) of Fig.4 and Figs. 6, 7 we plot the phase diagram obtained from minimizing Eq. (C8) for the band structure parameters of WS 2 , MoSe 2 and WSe 2 , respectively. Note that as opposed to the main text we do not specify the precise bilayer composition. For each TMD material here, one must consider a second "inactive" layer that induces the moiré potential and has band alignment properties that ensure it has a higher-in-energy conduction band. We use U = 60 meV and η = 1/2. Panel (a) and (b) are maps of the density of bands b and a, respectively, in the space of total density n and twist angle θ M . Panel (c) is the corresponding quasiparticle weight Z b . Panel (d) shows the relative filling at θ M = 3. There are two distinct regimes as a function of angle. For θ M < 3.5 the filling of band b is roughly split in half. Between n = 0 and n = 1 band b fills until it reaches a localized state (characterized by Z = 0 and n b = 1). Then band b fills continuously between n = 1 to 3. Finally, band b continues to fill until n = 4 is reached. The regime where band a is continuously filling realizes an orbitally selective Mott phase, where a Kondo lattice model is expected to be simulated with variable itinerant electron density.
On the other hand, for θ M > 3.5 the bands fill up oneby-one. In particular, band b fills completely between n = 2 and n = 4 with a Mott state at n = 5. Then above n = 5 it resets back into the Mott state and band a fills completely. This behavior thus, resembles a Stoner-like polarization of the species. However, we comment that the slave-rotor mean field tends to overestimate the size of band-polarized regions. This is because it overestimate the contribution of charge fluctuations to the interaction energy when the filing differs significantly from 1/2.

Two band slave rotor decomposition
As explained in the case of MoS 2 the electronic bands experience a much weaker spin-orbit coupling. Consequently, the shape and effective binding energies are almost identical [see Table I and Fig.3 (b)]. In this case it makes sense to decompose both bands in an unbiased manner where η = 1 and ∆¯ =¯ a −¯ b is much smaller than the band width and U .
Following the previous section, we now use four variational Hamiltonians of the form Eq. (7) and (C5), which generate a ground state |K τ , h τ , µ τ controlled by six variational parameters and subject to three constraints In panels (e)-(h) of Fig.4 of the main text we plot the phase diagram resulting from the minimization of the expectation value of Eq. (C9). As can be seen at n = 2 both Z a and Z b equal zero for small enough angles. The similarity of the bands of the two species leads us to propose electron doped MoS 2 as a candidate material to realize an SU(4) symmetric Mott insulator on a triangular lattice, which is an interesting problem on its own right which can is an interesting situation on its own right 54,55 .

Details of the numerical minimization procedure
In this appendix we provide the details of the numerical minimization procedure of Eq. (C8). To calculate this functional, we performed straightforward Brillouin zone integration on a square grid of size 150×150. The integration itself was performed by matlab's trapezoidal numerical integration, and the fermi-dirac distribution was written as: T being the inverse temperature. To broaden the discretization we use a finite temperature β = 60/max( a ). In addition, the minimal ground state that was found for the phase diagram in Fig.3 was found by matlab's minimization algorithm fmincon, which minimizes the functional Eq. (C8), under the constrains Eq. (C7) by means of the specified variational parameters. The optimization algorithm that was found to converge most efficiently was the interior-point algorithm, which is the default algorithm of fmincon. In the main text we have presented the results of a slave-rotor mean-field analysis, where bands of different species fill either one by one or simultaneously, depending on the twist angle. The latter scenario is is of particular interest to us as it gives way to the orbitally selective Mott phase.
Given that we have a number of unknown parameters, including η and the energy difference it is important to test the stability of the OSM state. In this appendix we compute a lower bound on the phase space volume of the OSM state in the space of ∆¯ and η.
To obtain this estimate we focus specifically on the commensurate filling n = 2 (or n = 4 for the strongly spin-orbit coupled bilaer WS 2 /WSe 2 ). If the bands fill one-by-one this filling point will be characterized by one completely filled band and another completely empty. On the other hand, if the bands fill simultaneously this filling value is likely to be characterized by a partial filling of both bands with total density of unity. Here we will make a restrictive assumption, that both bands are at filling unity, where one is in a Mott state and the other is either metallic or also in a Mott state. The phase diagrams we have computed are consistent with this behavior at filling n = 2 (or 4).
Under this assumption we can estimate the expectation value of the Hamiltonian Eq.6 within these restrictive trial states: and compare which of them has a lower energy. The first two states represent fully polarized states and thus corresponds to the scenario where the band fill one-byone. The third state however, is where the filling is shared between the two bands. We will further assume that band a is the flatter of the two bands, and is in a Mott state at n a = 1.
The energy per unit cell of these three states is given by Here Ω is the total number of sites. T it ≤ 0 is the sum of the negative kinetic energy associated with half-filling band b, and the interaction energy associated with the charge fluctuations ∆n 2 b = (n b − 1) 2 at the half-filling point. Thus, when T it < 0 band b remains in a metallic state and reaches zero when it falls into a Mott state as well.
When E ab < min(E a , E b ) the third state (simultaneous filling) is more favourable energetically. On the other hand, when E a or E b are minimal, the ground state can be band-polarized (but not necessarily), where the bands fill one-by-one.
Comparing these energies, we conclude that the partial filling state is stable at least in the regime In Fig. 8 we plot the stability region defined by Eq. (D4) for the case of T it = −U/4. The case of T it = 0 is also plotted for comparison (dashed line). Given that T it ≤ 0, there exists such a regime for any value of η. Note that the width in ∆¯ of this window scales with U at strong coupling. We thus conclude that the regime of partial occupation of both bands is wide and robust to parameters such as the ratio of species interaction η and splitting of the bands ∆¯ .
In Fig. 9 we plot the width of the OSM phase ∆n = n max − n min averaged over the angles 2 • and 5 • , which is obtained from the numerical minimization of Eq. (C8). Here n max and n min mark the boundaries of the OSM phase per angle (maximal value is 2), as shown in Fig.  4 .(c). The white dashed line is the analytic estimate Eq. (D4). We thus, conclude that the OSM phase is not sensitive to parameters and is a generic feature of the phase diagram of electron doped moiré TMDs.

E. Spin Exchange interactions and expected phenomenology
Inside the OSM phase charge fluctuations of the localized species are quenched and therefore, the relevant inter-species interactions are of spin-exchnage type. In this section we discuss two such interactions and their possible influence on the magnetic ground state.
We anticipate that the ferromagnetic Hund's coupling is the largest exchnage mechanism where S i is the spin of the localized moments at site i (let us assume they belong to species a) . Using standard harmonic oscillator states, we estimate J H ≈ 0.2U . Upon approaching the meting point of the OSM state however, other interaction terms that compete with Eq. (E1). For example, the interaction that scatters across the original Brillouin zone. To see this let us first write this term in basis of the original operators Eq. (1) H J = J P k,k ,p c † k+p a↑ c † k −p a↓ c k b↓ c k b↑ + h.c., where J P /U ∼ a/a M ≈ δ 2 + θ 2 M (for an angle of θ M = 3.5 • we obtain J P ≈ 0.1U ). Taking into account the moiré potential and the slave-rotor decomposition described above, this interaction assumes the form H J =J P k,k ,p f † k+p a↑ f † k −p a↓ ψ k b↓ ψ k b↑ + h.c. (E2) whereJ P = Z a J. Thus, when the band a is localized and Z a = 0 this term vanishes. However, if the local moments are incorporated into the Fermi surface (through the formation of a heavy Fermi liquid) they reacquire a finite quasi-particle weight Z 49 . Thus, this interaction can become important close to the melting point of the local moment lattice. As mentioned above, inside the OSM phase we expect the dominant exhnage to be Hind's and therefore spincorrelations to be ferromagnetic as in Ref. 30. When the Hund's coupling Eq. (E1) is dominant the main influence of the itinerant electrons inside the OSM state is to mediate long-ranged RKKY interactions 58 where ν 0 is the density of states. Thus, as the filling of the itinerant band is modified from zero to 2, the nearestneighbour interaction mediated by the electrons can be tuned from ferromagnetic in the dilute limit to antiferromagnetic and back to ferromagnetic (going through a van Hove singularity). This interaction is added to the direct superxchnage between sites 19 , which may have a cooperative effect or frustrate the magnetic interactions.

F. The heavy Fermi liquid state
When the antiferromagnetic correlations dominate we expect a heavy-Fermi liquid state to compete with the internal magnetic interactions. For completeness, in this appendix we compute the Kondo temperature assuming J P = 4 meV within the large-N mean-field theory 65 . We find hat T K ∼5-15 K inside the OSM phase and therefore, we expect that if the antiferromagnetic correlations dominate the OSM melting transition can be accompanied by a detectable heavy Fermi liquid state. Moreover, in this case the tunability of the itinerant electron density may allow to tune though the Doniach phase diagram 26 .
Let us breifly describe the large-N mean field theory. The dispersion of the two species is taken to be where the density of each species n τ is set separately using the Lagrange multipliers µ τ according to their values in the slave-rotor mean-field calculation. We then decouple Eq. (E2) using the mean-field hybridization χ = J/2 f † bσ ψ aσ + c.c.. We then solve for χ self-consistently, while tuning µ a and µ b to conserve the density of a and b on average. The Kondo temperature is then estimated by seeking the lowest temperature where the self-consistent solution for χ = 0. To obtain the Kondo temperature T K we estimate the lowest temperature where χ = 0.
In Fig. 10 we plot the resulting Kondo temperature as a function of twist-angle and density forJ = 4 meV. As can be seen, the Kondo temperature measurable in standard cryosthetics and might be physically important.