Doping a Wigner-Mott insulator: Exotic charge orders in transition-metal dichalcogenide moir\'e heterobilayers

The moir\'e pattern induced by lattice mismatch in transition-metal dichalcogenide heterobilayers causes the formation of flat bands, where interactions dominate the kinetic energy. At fractional fillings of the flat valence band, the long-range electron interactions then induce Wigner-Mott crystals. In this Letter we investigate the nontrivial electronic phases appearing away from commensurate fillings. Here, competing phases arise that are either characterized as doped Wigner-Mott charge transfer insulators or alternatively, a novel state with frozen charge order yet is conducting: the 'electron slush'. We propose that an extremely spatially inhomogeneous local density of states can serve as a key signature of the electron slush.

When two atomically thin materials are intentionally misaligned or have a lattice mismatch, a long-range geometric moiré pattern emerges.This results in a drastic reduction of the electronic kinetic energy, paving the way for new strongly correlated phases.[1,2] Though the first signs of correlations were seen in twisted bilayer graphene [3,4], it is particularly appealing to use instead monolayer semiconductor transition metal dichalcogenides (TMDs) MX 2 where M = W, Mo and X = S, Se, Te. [5][6][7][8][9][10] In TMD moiré bilayers, the effective interaction strength U scales as the inverse moiré length U ∼ a −1 M whereas the flat band kinetic energy at the top of the valence band scales as W ∼ a −2 M .Consequently, this back-of-the-envelope argument suggests U/W ∼ a M , meaning there is no limit as to how strongly coupled a TMD bilayers can be (Appendix A)! Furthermore, combining various TMD monolayers with or without a relative twist angle allows for a wide degree of variability in constructing the physics of flat bands [11].Most notably, in TMD hetero-bilayers (combining two different TMDs) the electronic states are localized in only one layer and form an effective triangular lattice Hubbard model [12][13][14][15][16][17].
Exactly at these fractional fillings, the Wigner-Mott crystal is an insulator, which can be described on the mean field level [13,[33][34][35][36][37].Away from fractional fillings, the situation becomes less clear.While there is a large theoretical literature on interacting triangular lattice models, most focus on superconducting instabilities [38][39][40][41][42], which is not observed in experimental moiré systems.Instead, a realistic possibility is the emergence of compet- 1.The local density of states (LDOS) of a typical "electron slush" phase at incommensurate fillings of a TMD heterobilayer, obtained using non-charge-self-consistent Hartree+CPA+DMFT, with T = 1.2K, n = 1/2, and U/t = 38.3.Here, the electrons are frozen in an amorphous configuration.The presence of screened Coulomb interactions causes a soft spectral gap (Fig. 4) and nonzero conduction (Fig. 3).However, the LDOS is spatially extremely inhomogeneous, and dominated by only a few peaks at the Fermi level (a).At lower energies, here at E = −20 meV (b), the LDOS shows more structure of the amorphous slush state, which can be observed using scanning tunneling microscopy.
ing charge ordered phases, as was recently suggested by Monte Carlo simulations [43].Alternatively, the system can become a doped Wigner-Mott insulator [44] where the charge order is given by the nearest fractional filling.
Aside from variational approaches close to integer filling [45], studies of the interplay of Mottness and charge order beyond mean-field theory are all based on versions of DMFT, augmented with either self-consistent Hartree-Fock for the charge order [21-23, 46, 47] or GW for the momentum-dependent self-energy.[48,49] Following the former approach, at most incommensurate fillings we identify a novel metastable "electron slush" phase.Here the electron charge freezes spontaneously in amorphous real-space patterns, whose signatures in the local density of states are shown in Fig. 1.Despite the frozen charge order, the screening of long-range interactions prevents the opening of a spectral "Coulomb" gap.In addition, the random electrostatic fields produced by such amorphous charge order induce local doping of the Wigner-Mott solid, resembling the effect of site disorder in ordinary Mott insulators.As a result, the system remains a conductor -albeit a very bad one, with an extremely high resistivity surpassing the Mott-Ioffe-Regel limit.This stands in great contrast to competing "doped Mott-Wigner insulator" phases, where the fractional periodic charge order remains, and the dopants realize a strongly renormalized (but highly conducting!) Fermi liquid at low temperatures.Physically, this amorphous charge order can be seen as a "compromise" between several closely competing periodic Wigner orders, at incommensurate fillings.The strange spectral and transport properties characterizing such amorphous Mott matter are the central discovery of this work.We emphasize that the frozen charge order is the result of very strong nonlocal Coulomb repulsion V /t ≫ 1, which can be captured self-consistently on a semi-classical Hartree level.On top of that, we use DMFT to gain insights into the effect of such charge freezing upon the spin and transport effects associated with strong electronic correlations due to the Hubbard U .
Moiré model -To model the aligned WS 2 /WSe 2 hetero-bilayers we follow the continuum formulation proposed by Wu et al. [12].We briefly summarize their approach here, a more detailed description is provided in Appendix A. The continuum model starts with the top of the valence band in monolayer WSe 2 , where the holes have an effective mass of m * = 0.36m e .
Aligned (θ = 0) WS 2 /WSe 2 heterobilayers have a moiré length of approximately a M = 7.9 nm.To obtain the strength and phase of the periodic moiré potential, we performed density functional theory calculations using the approach of Refs.[12,50,51].For WS 2 /WSe 2 bilayers we obtained (V, ψ) = (7.7 meV, −106 • ).Note however, that the accuracy of the magnitude V is debatable.Recently, STM measurements revealed a moiré potential in excess of 100 meV [52], whereas theoretical work suggests a complete absence of a moiré potential [16].Additionally, the magnitude is very sensitive to the precise interlayer distance (and can thus be tuned using pressure as was done in twisted bilayer graphene [53]).Therefore, one should bear in mind that the provided values for the electronic bandwidth can be significantly different in real systems.
Nevertheless, with this value of moiré potential a flat band emerges that can be represented by the extended Hubbard model on a triangular lattice.The nearest neighbor hopping is t = 1.9 meV, and we ignore longerranged hoppings.The corresponding Wannier orbitals [54] are approximately Gaussian around the W/Se stacking center with width σ = 0.125a M .We use the screened Coulomb potential to calculate the effective interaction strengths.With a screening length of d = 20 nm, this provides an onsite repulsion of ϵU = 1.7 eV, and a nearest neighbor repulsion of ϵV = 0.6 eV.Even with quite large values for the dielectric constant the onsite repulsion is stronger than the bandwidth, and indeed even V ≫ t!Because both the bandwidth (via the moiré potential and interlayer coupling) and the interaction strength (via engineering of the dielectric environment) are tunable, in the remainder of this paper we vary the ratio U/t.Our choices are U/t = 38.3,19.1 and 12.8, corresponding roughly to dielectric constants of ϵ = 23.4,46.8 and 70.2, respectively.For the strongest choice of U/t, the corre-sponding nearest neighbor repulsion becomes V /t = 13large enough to induce nontrivial charge order.
Charge order -Having established that the long-range repulsion is quite strong, we next discuss the development of charge order.In the absence of a lattice potential, a dilute two-dimensional electron gas will form a triangular Wigner crystal [55,56].Whenever an underlying triangular lattice is present, the triangular Wigner crystal can only be formed at fractional fillings given by n = with m i integers.The largest fractions are n = 1/3, 1/4, 1/7, and 1/9; one can also have a generalized Wigner crystal of holes at fillings n = 2/3, 3/4, etc.At filling n = 1/2 the long-range repulsion favors stripe order [57].At all of these fractional fillings a correlated insulator state has been observed in WSe 2 /WS 2 [24][25][26][27].
Because the onsite repulsion U is large, to leading order we can project out double occupancy.As a first step to identify various forms of charge order across the phase diagram, we use self-consistent Hartree theory for long-range interacting spinless electrons, with screened Coulomb interaction V (r) = V0 r e (r−1)/d where r is measured in moiré lattice constants, d = 20 nm the screening length and we choose V 0 = 30 meV (Appendix B).The zero-temperature free energy of various charge ordering patterns is shown in Fig. 2, left.Here we compare the free energy of five different possible charge patterns as a function of filling n.The most robust of these configurations is the Wigner-Mott crystal with 3 sublattice sites, yielding a triangular Wigner crystal at n = 1/3 and a honeycomb crystal at n = 2/3.Similarly structures appear at n = 1/4, 1/7, and other fractional fillings.In addition, we confirm the presence of the n = 1/2 stripe phase [27,57].
Throughout the entire range of filling 0 < n < 1, we further compare these Wigner crystal states with aperiodic frozen charge patterns, also known as amorphous states.For this we used a large 12×12 unit cell, equivalent to 144 moiré unit cells, and found a local minimum of the Hartree free energy corresponding to an inhomogeneous charge distribution.While most electron crystal states are stable upon small doping, there appears in between each rational fillings a regime where amorphous patterns are the most stable ones.If the system is cooled fast enough and with infinitesimal disorder, however, crystalline orders or macroscopic phase separation can be avoided at all fillings.The electrons then freeze into an amorphous solid, behavior that has been observed in θ-organic materials [58][59][60], consistent with recent theory [57].The emergence of similar metastable amorphous structures has also been reported in various TMD layered materials featuring Wigner crystals and glasses made of "star of David" polarons [61,62].
The results presented in Fig. 2 were obtained for the interaction strength V /W ≈ 2 -nearest neighbor repulsion is about twice the noninteracting bandwidth.This ratio is highly dependent on the dielectric environment and the strength of the moiré potential.With the parameters used here we find that the n = 1/3 Wigner crystal is stable up to T c ≈ 12 meV, which is a factor 2-3 higher than observed in experiments [24,26].It is therefore likely that in real materials V /t ≈ 5, which would correspond to U/t ≈ 15.Since mean field theory systematically overestimates the tendency to order, these approximations for U/t and V /t are lower bounds.
Doped Wigner-Mott insulators -Having established various forms of charge order present, we now investigate their impact on Mott-Hubbard correlations, within a setup based on dynamical mean-field theory (DMFT) [63,64]; a detailed description of our calculation methods can be found in Appendix C. Representative results for the resistivity as a function of temperature are shown in Fig. 3, and the corresponding spectral functions are shown in Fig. 4.
Our starting point is the Wigner-Mott insulator, such as the state we find at n = 1/3 filling.In our system, the on-site repulsion U is generally larger than the gap induced by the charge order, so the resulting Wigner-Mott state can therefore be recognized as a charge-transfer insulator [13,65,66].Here the electronic states below the Fermi level are Mott-localized on the occupied sublattice, whereas the states above the Fermi level reside on the unoccupied honeycomb-like sublattice.This is illustrated in Fig. 4b, where for n = 1/3 we indicate the charge transfer (CT) gap, and the upper and lower Hubbard band (UHB/LHB).
Exactly at fractional filling, the strong onsite repulsion leads to three distinct transport regimes, shown in Fig. 3b.At high temperatures above the crystalline ordering transition T > T c we find a typical ρ = A + BT linear resistivity bad metal regime with resistivity comparable or even above the Mott-Ioffe-Regel limit.Just below the charge-ordering transition (T ≤ T c ) the resistivity jumps up and the system becomes insulating.When the interactions are weaker, in our model for U/t = 12.8, we find an intermediate situation where the charge order can form but the onsite U is not strong enough to open a Mott gap, leading to low T Fermi-liquid behavior with well-defined quasiparticles.Note that this is in contrast to Ref. [37] where it is assumed that the Mott and the charge order transition always happen simultaneously.In our picture, as in previous work on Wigner-Mott systems [21][22][23], the weakly first-order Mott-like metal-insulator transition occurs within the charge-order state.Our scenario is actually realized in many real materials [61,62] featuring Wigner crystallization in lattice systems, which can lead to either insulating or metallic charge-ordered ground states.
Let us now turn our attention to doping the Wigner-Mott state.For weaker interactions (U/t = 12.8), doping a metallic charge density wave retains the conventional Fermi liquid behavior.However, when the onsite Mott correlations are large we see different behavior, depending on whether we dope with holes or electrons.
For electron-doping (n = 1/3 + δ with δ > 0, Figs.3c  and 4c), the dopants appear in an otherwise empty honeycomb sublattice.As a result, the dopants are moderately interacting and form pinball liquid [46,47] with a renormalized Fermi liquid regime at low temperature.We find that the associated Fermi liquid coherence temperature vanishes linearly as we approach Wigner-Mott insulator, T F ∼ |n − n c |.This is a clear prediction of our theory that can be tested in experiment, similar to recent measurements of bandwidth-controlled Mott criticality [19,67].
This Fermi liquid regime evolves at intermediate temperature into an incoherent high-resistivity phase.Here, the shifting of the chemical potential into the charge ordering gap competes with the loss of quasiparticle coherence, resulting in non-monotonic behavior of the resistivity.Finally, as is the case at fractional filling, when charge order is lost the high-temperatures, linear-T bad a. n=0.317  3, is due to a shifting chemical potential which can be clearly seen in the left figure.Upon increasing temperature the charge-transfer gap closes and a strange metal phase appears.The DOS for the amorphous state (d) has similar features (LHB, CT gap, UHB) but its features are smeared out due to the inhomogeneous charge patterns, leading to a soft gap in the spectrum and consequently bad metal behavior (Fig. 3).
metal behavior is recovered.
Upon hole-doping (n = 1/3 − δ, Figs.3a and 4a), the dopants are populating the sites of the triangular Wigner-Mott lattice.While at extremely low temperatures a FL regime does appears (not shown), for most temperatures the Fermi level lies within the incoherent lower Hubbard band, leading to bad metal behavior.As the temperature increases, the chemical potential shifts into the charge-transfer gap leading to alternating weakly insulating and weakly metallic states.More details of the transport behavior within the entire regime of dopings and temperatures around the Wigner-Mott regime are shown in Appendix D.
Electron slush -As shown in Fig. 2, however, for the largest part of the phase diagram, at most incommensurate fillings, the system is likely to freeze in an amorphous charge pattern.Because any amorphous charge order results in a random distribution of internal electrostatic potentials, within our DMFT setup the rest of the calculation reduces to solving an appropriate Hubbard-like (i.e.charge-transfer) model supplemented with random site energies, which we tackle within the well-known CPA-DMFT approach [68] (see Appendix C and E).In the following we focus on the amorphous states at n = 1/2 and for U/t = 19.1,shown in Figs.3d and 4d; within our setup, the qualitative behavior at other fillings and other (comparable) values of U/t is essentially the same.Note that n = 1/2 is a fractional and not an incommensurate filling; we performed our DMFT calculations at this filling for convenience; the properties of the electron slush are the same at fractional and incommensurate fillings.
Even though the electrons clearly develop local frozen charge order, this does not lead to the opening of a hard gap in the spectrum.This situation reminds us of the formation of a soft "Coulomb" gap in electron glasses [69][70][71][72].However, here the screening of the long-range Coulomb interaction (due to the presence of gate on the TMD heterobilayer) causes a further "filling" of the soft gap.As a result, despite the frozen charge order, there is no gap towards transport: the amorphous ordered state can still conduct, although they do so very poorly.Our DMFT results indicate at low temperatures the formation of weakly-developed quasiparticles around the Fermi level, with spectral features reflecting the amor-phous charge background in question.This combination of frozen charge order with motion of electrons reminds us of the popular "slush" drinks which are liquids (conducting) yet frozen.We therefore propose to call this novel phase an electron slush.
Note that the amorphous charge order of the electron slush can be characterized by the presence of shortrange stripe correlations, which has also been observed in Monte Carlo simulations [26,43,57].Therefore, scanning tunneling microscopy (STM) [28,29] topography should reveal the frozen charge configurations.A measurement of the local density of states (LDOS), however, leads to an interesting local realization of the soft gap feature.As shown in Fig. 1, the LDOS at the Fermi level is dominated by a few localized peaks at large distances, whereas at energies away from the Fermi level the amorphous structure is more spread out.This stark difference between topography and LDOS is another predicted feature of the electron slush phase, which awaits experimental verification.
Outlook -Summarizing, we show that, for most fillings, an amorphous charge ordered metallic state -the electron slush -appears in strongly correlated TMD heterobilayers.Our results suggest to perform transport and STM experiments on TMD heterobilayers to study the interplay between doping a Wigner-Mott insulator and amorphous charge order.
Note that in this Letter we ignored the role of spinorbit coupling (SOC) in the moiré flat bands [17].Also, in this work we did not include the role of quenched disorder.It is known that TMD mono-and bilayers have significantly more disorder than other Van der Waals materials such as graphene, notably in the form of vacancies.[73] Since disorder strongly affects the zerotemperature limit of the resistivity, our transport results only apply at intermediate to high temperature.On the other hand, disorder is likely to further stabilize the electron slush phase since it promotes amorphous charge configurations.Our work shed light on the interesting but mysterious amorphous phase, which could very well play a central role in strongly correlated TMD heterobilayers at low band filling.Here we took the first steps to theoretically describe the likely features of this exotic regime, although more accurate treatments of the interplay between strong correlations and induced/effective disorder effects may be necessary, especially at the lowest temperatures.This research direction remains a challenge for future work.
A particularly insightful way to understand the flat bands in heterobilayers is the continuum model proposed by Wu et al. [12].Here we apply their method to WS 2 /WSe 2 heterobilayers.
In a single TMD monolayer the top of the valence band has a parabolic dispersion around the K and K ′ point.Strong spin-orbit coupling causes a spin-valley splitting, so that the states at K/K ′ are only singly degenerate and carry opposite spin.The band-structure in a single valley can thus be approximated as where q = k−K/K ′ and the effective mass m * = 0.36m 0 for WSe 2 .
Let us now look at a bilayer system.A monolayer WSe 2 has a lattice constant of approximately a 1 = 3.325 Å, whereas monolayer WS 2 has a 2 = 3.191 Å. Therefore even when the two layers are aligned, a moiré pattern emerges due to the lattice mismatch.The length scale of this pattern is given by where θ is a twist angle.When the twist angle is small, and writing a 2 ≡ a 1 (1 − δ), we obtain An aligned (θ = 0) WS 2 /WSe 2 heterobilayer therefore has a moiré length of approximately a M = 7.9 nm.Note that only heterobilayers where the chalcogenides are the different in the two layers lead to a Moiré pattern when aligned.Furthermore, the Moiré length is maximal for our combination of S and Se: WS 2 /WTe 2 has a M = 3.2 nm whereas WSe 2 /WTe 2 has a M = 5.0 nm.This makes the platform WS 2 /WSe 2 likely the most correlated among all aligned TMD heterobilayers.
Electronically, monolayer WS 2 has a larger bandgap than monolayer WSe 2 , and the bilayer will have a type II band alignment of the two layers.This means that the top of the valence band of the bilayer consists of electronic states confined to the WSe 2 layer only.Note that a perpendicular electric field can change the band alignment of the two layers, thus bringing the WS 2 valence band at the same energy as the WSe 2 valence band, as was done recently in MoTe 2 /WSe 2 bilayers.[20] Here, however, we study the WS 2 /WSe 2 bilayer in the absence of an electric field.In that case, the electronic states in the WSe 2 feel a position-dependent Moiré potential ∆(r) due to the presence of the WS 2 layer.
The Moiré potential can be approximated using plane waves where G M i with i = 1 . . .6 are the six reciprocal Moiré vectors.Without loss of generality, we set G M 1 = 4π 3a M (1, 0) and the other five are just sixfold rotations of the first reciprocal vector.Because ∆(r) must be real, and it has three-fold rotational symmetry, we have [12,50,51] which means we can parametrize the Moiré potential using only two parameters (V, ψ) The actual magnitude of the Moiré potential can be estimated using ab initio density functional theory.For this we use Quantum Espresso [74,75] with a Coulomb cut-off [76] to reproduce the two-dimensional nature of the heterostructures.We used a lattice structure where the interlayer distance, as measured by the z-distance between the W atoms, is d W −W = 6.6 Å, and the z-distance between the W and chalcogenides is 1.65 Å.The idea, following Refs.[12,50,51], is to calculate the energy of the top of the valence band in a small unit cell WS 2 /WSe 2 bilayer where the top layer is shifted with a displacement d.In the full Moiré unit cell, the Moiré potential follows the same energy dependence as the top of the valence band in the small unit cell.
For aligned WS 2 /WSe 2 bilayers we obtained (V, ψ) = (7.7 meV, −106 • ), see Fig. 5a.This corresponds to a Moiré potential peaked where the Se atoms are above the W atoms in the other layer, consistent with scanning tunnelling microscopy (STM) measurements.[15] Note however, that the accuracy of the magnitude V is debatable.Recently, STM measurements revealed a Moiré poten-tial in excess of 100 meV, [52] whereas theoretical work suggests a complete absence of a Moiré potential.[16] Additionally, the magnitude is very sensitive to the precise interlayer distance (and can thus be tuned using pressure as was done in twisted bilayer graphene [53]).Therefore, throughout this work, one should bear in mind that the provided values for the electronic bandwidth can be significantly different in real systems.
The potential Eq. ( 4) causes backfolding of the bandstructure Eq. ( 1) into the mini-Brillouin zone.At the edge of the mini-Brillouin zone the potential opens up a gap, leading to the formation of flat bands.The bandstructure of aligned WS 2 /WSe 2 is shown in Fig. 5b.Since we are dealing with a single band (per valley), a straightforward Fourier transform provides us with the tight-binding parameters.The nearest neighbor hopping is t = 1.9 meV, and the next-nearest neighbor hopping t ′ = −0.4meV.Longer-ranged hopping parameters fall off exponentially and can be safely ignored.
Similarly, Wannierization [54] using a projection onto a Gaussian wavepacket ansatz yields the Wannier orbitals which, in this case, reduces to approximately Gaussian around the W/Se stacking center with width σ = 0.125a M .This simple shape of the Wannier orbital allows us to directly extract the onsite repulsion U , and longer ranged Coulomb interaction parameters V .For this, we use the standard screened Coulomb potential from having a screening layer at distance d from the TMD heterobilayer.With a screening length of d = 20 nm, this provides an onsite repulsion of ϵU = 1.7 eV and a nearest neighbor repulsion of ϵV = 0.6 eV.Even with quite large values for the dielectric constant the onsite repulsion is stronger than the bandwidth, and indeed even V ≫ t.Because both the bandwidth (via the Moiré potential and interlayer coupling) and the interaction strength (via engineering of the dielectric environment) are tunable, in the remainder of this paper we vary the ratio U/t.Our choices are For the strongest choice of U/t, the corresponding nearest neighbor repulsion becomes V /t = 13 -large enough to induce nontrivial charge order.
Note that the size of the Wannier orbital changes upon changing the moiré length.In the literature, both σ ∼ a M [77] and σ ∼ √ a M [12] have been reported.The latter limit is valid in the limit where the moiré potential V is larger than the kinetic energy associated with the reciprocal moiré lattice vector ℏ 2 G 2 2m * , whereas the linear dependence is valid when V is small.For untwisted WS 2 /WSe 2 , we find ℏ 2 G 2 2m * ≈ 87 meV whereas V = 7.7 meV, clearly placing WS 2 /WSe 2 in the regime where the Wannier orbital size is proportional to the moiré length, σ ∝ a M .

APPENDIX B: SPINLESS HARTREE METHODS
The Hartree self-energy is defined as Σ H (r) = r ′ V (r − r ′ )⟨n r ⟩ where the Coulomb interaction is approximated as V (r) = V0 r e (r−1)/d where r is measured in Moiré lattice constants, d = 20 nm the screening length and we choose V 0 = 30 meV.Including the Hartree selfenergy, the total free energy is given by which is iteratively minimized, starting from a random initial configuration at low T = 0.1 meV, on various large supercells containing multiple moiré unit cells.
To obtain the free energy of commensurate charge ordering patterns (Fig. 2 of the main manuscript), we take small supercells containing 2, 3, 4 or 7 moiré unit cells.For these supercells there is a unique configuration that yields the lowest free energy.For example, for the curve corresponding to "1/3" order, we calculated the Hartree solution where the lattice symmetry is broken to a larger supercell containing 3 moiré unit cells.This means the charge density can be different on the three moiré unit cells contained in the supercell; let's call them n 1 , n 2 , and n 3 .The average filling is n = (n 1 + n 2 + n 3 )/3.At exactly n = 1/3 filling, the solution exists where n 1 = 1 and n 2 = n 3 = 0, which is shown in the side-panel of Fig. 2 of the main manuscript.At other fillings, the "1/3 order" curve corresponds to the lowest free energy possible given any choice of n 1 , n 2 , and n 3 .
To obtain the free energies at incommensurate charge orderings, we considered supercells containing L × L moiré unit cells with L = 8, 12, 24 and 36.For each system size we started with a random initial configuration and found a (meta)stable mean field solution at zero T .After that, we increased temperature starting in each configuration to find the evolution of the charge order at finite temperature.
There are exponentially many different metastable configurations that can be found in this way.An indication of this is shown in Fig. 6 for L = 12.The onsite Hartree energies act as input for the subsequent DMFT calculation, for which we calculated 100 metastable amorphous configurations at L = 12.

APPENDIX C: DYNAMICAL MEAN FIELD THEORY METHODS
Following the previous section, we know that the flat valence band in TMD Moiré materials can be treated as an extended Hubbard model on a triangular lattice where c † i,σ , c i,σ are the fermionic creation/annihilation operators, i, j are site labels to be summed over, n iσ ≡ c † iσ c iσ the density operator, t ij the hopping amplitude, V ij the intersite long range Coulomb interaction strength and U the onsite Coulomb interaction strength.The hamiltonian Eq.( 6) differ from the single-band Hubbard model by considering also intersite correlations, which can be reduced to a local form using Hartree theory by replacing the density operator n j,σ with the spinless expectation value ⟨n j ⟩ i,j̸ =i,σ1,σ2 where the effective site energy ε i is written as By using the Hartree approximation Eq. ( 7), all interactions in the extended Hubbard model Eq. ( 6) are now local In the framework of single-site DMFT, Eq. ( 9) is solved first by identifying the corresponding Anderson impurity problems then subject to DMFT self-consistency conditions on the bath ∆ (ω).In below we will outline the general single-site DMFT procedure to solve Eq. ( 9) for all commensurate fillings.
Except at n = 1, the corresponding Anderson impurity problem consists of m > 1 individual sites in the unit cell, depending on the Wigner crystal configuration in Fig. 2 and hence we have ∆ i=1,2,...,m (ω) → AIM solver → Σ i=1,2,...,m (ω) , (10) which is then subject to m DMFT self-consistency equations Specifically, in Eq. ( 10), we use real frequency Iterative Perturbation Theory (IPT) for arbitrary filling as the Anderson Impurity Solver(AIM solver) [78].The term G ii (ω) is computed from the local Green's function G (ω) where η is the broading term, I is the identity matrix, E (k) is the dispersion matrix and Σ (ω) is the diagonal selfenergy matrix A solution is obtained iteratively by: 1. starting from initial ansatzes ∆ i (ω), 2. solving the corresponding Anderson impurity problems to obtain Σ i , 3. use selfconsistency condition Eq. ( 11) to calculate ∆ i (ω) for the next iteration.Specifically for the triangular lattice at 1/3 filling in this paper, the exact form of the band structure is the following: where ab + e ikδ (2)  ab + e ikδ (3)   ab ; ac + e ikδ (2)  ac + e ikδ (3)   ac ; bc + e ikδ (2)  bc + e ikδ (3)   bc ; and and a is the lattice spacing.
In the case of an amorphous moire lattice, there is no periodicity hence the system cannot be reduced to repeating unit cells (for n = 1/2 see Fig. 2).Calculations indicate that the Hartree site-energy is essentially randomly distributed and follow a certain probability distribution P (ε), shown in Fig. 7.A possible mean-field description of such a system is the Coherent Potential Approximation (CPA).In CPA, spatial variation is disregarded, such that the impurity within the DMFT framework is replaced by the average impurity, described by average Green's function Ḡ (ω) where − 1 π ℑG (ϵ, ω) gives the LDOS depending on the local Hartree site-energy.Eq. ( 18) leads to the effective self-energy The DMFT self-consistency equation is closed by where and E k = −2t cos (k x ) + 2cos √ 3/2k y cos (1/2k y ) .To sum up, the DMFT+CPA recipe for amorphous system is as follows: 1. Input hybridization function ∆ (ω); 2. Solve Σ (ω, ε) for each Hartree site energy ε, using AIM impurity solver (Eq.( 10)); 3. Compute average Green's function Ḡ (ω) (Eq.( 18)); 4. Compute effective self-energy Σ (ω) (Eq.( 19)); 5. Find ∆ (ω) for the next iteration and repeat from 1 till converge (Eq.( 20)).
To get a good statistic, in step 2, we need to sample at least 100 Hartree site energies (solve for 100 impurity problems).Each DMFT loop takes 30 iterations on average, and on top of the DMFT loop, in step 6, the chemical potential µ needs to be adjusted 5 ∼ 10 times to keep n fixed.So, to sum up, for one histogram P (ϵ), we need to solve roughly 100 * 30 * 5 = 15K impurity problems.
To obtain the LDOS mapping of large area shown in Fig. 1 and Fig. 8, one firstly obtains how exactly the electrons freeze in an amorphous charge ordered pattern, in other words, the exact location (x i , y i ) of each Hartree site-energy ϵ i , via minimizing the total free energy on multiple Moiré unit cells.Then the LDOS for the states at energy ω is where δ 2 x , δ 2 y are chosen to be 0.12a M , so that the width of the Gaussians is approximately the width of the Wannier functions, which is 0.125a M .For transport calculations, the DC conductivity is calculated via Kubo formula [79].: with We only focus on the xx component of the conductivity.For amorphous state, E (k) is replaced by ϵ (k).
In our results (Fig. 9), we find that in the homogeneous case a Wigner-Mott insulating state around the commensurate filling at low temperature, while at high temperature charge-order is lost and the system becomes that of a bad metal with resistivity exceeding that of the Mott-Ioffe-Regel limit.We also investigate the effects of doping away from the commensurate filling, which results in a fermi-liquid at low temperatures and bad metal at higher temperatures.Exactly at n = 1/3 filling, the system becomes insulating due to a combination of charge order and Mott localization leading to a large charge-transfer gap.At temperatures where the charge order vanishes, linear resistivity is found indicating bad metal behavior.Electron or hole doping yield different phases due to the asymmetric nature of the charge-transfer gap.Upon electron doping, the low temperature phase is a heavy Fermi liquid with T F L vanishing linearly with doping.On the hole doped side, the Fermi liquid is pushed to extremely low temperatures.f = 1 state.We also extract the activation gap from the Arrhenius plot.We can see the activation gap smoothly goes down but it extrapolate to zero at lower U/t than where the quasi-particles arise (Red star in Fig. 10(c)), since it's a weakly first order transition.

APPENDIX E: STABILITY OF THE ELECTRON SLUSH
In this section we discuss whether the electron slush is a stable phase if we would include higher order corrections to our approximation scheme for the electron self-energy.Note that in the classical limit of very small t, the existence of amorphous metastable states is well established.This has been discussed in, for example, Ref. [72].Such amorphous metastable states on a triangular lattice with Coulomb interactions have even been observed in organic compounds, see Refs.[59,60].
This charge order is stable with respect to the introduction of a quantum hopping t.As such, introducing higher order corrections to the electron self-energy will inevitably lead to a reduction of the charge order but not to a complete suppression.Since charge order has been very clearly observed in WS 2 /WSe 2 , see Refs.[24][25][26][27][28][29][30][31][32], we are certainly in a range where nonlocal corrections to the electron self-energy affect quantitatively but not qualitatively the charge order.
What is left is the interplay between the developed charge order and the local onsite correlations as captured by DMFT.The question is whether the Hartree static mean field theory for the charges is affected by the nonstatic local self-energy.
To check this, we calculated the occupation number n versus different Hartree site energies ϵ for the amorphous state for T = 1.2 K, see Fig. 11, and compare them with the occupation numbers after taking into account the DMFT solutions.The Hartree site energies were extracted from the distribution P (ϵ) as shown in Fig. 7.Note that the average occupation number equals n = dϵ P (ϵ) n(ϵ) is 1/4 in the spinfull model.
In Fig. 11 we see that while the DMFT rounds the charge occupations in a ∼ 10 meV window, it does not qualitatively change the unequal site occupations characteristic of amorphous charge order.The randomess of n Hartree mainly comes from the randomness of a specific bath chosen in each calculation.This is consistent with earlier works that fully self-consistently contained both onsite U and nearest-neihbor V on a square lattice [21]: the inclusion of large U does not change the charge order transition.
Nevertheless, more quantitative accuracy is expected when one uses a technique that also includes the nonlocal dynamic self-energies, such as EDMFT.
FIG.1.The local density of states (LDOS) of a typical "electron slush" phase at incommensurate fillings of a TMD heterobilayer, obtained using non-charge-self-consistent Hartree+CPA+DMFT, with T = 1.2K, n = 1/2, and U/t = 38.3.Here, the electrons are frozen in an amorphous configuration.The presence of screened Coulomb interactions causes a soft spectral gap (Fig.4) and nonzero conduction (Fig.3).However, the LDOS is spatially extremely inhomogeneous, and dominated by only a few peaks at the Fermi level (a).At lower energies, here at E = −20 meV (b), the LDOS shows more structure of the amorphous slush state, which can be observed using scanning tunneling microscopy.

FIG. 2 .
FIG.2.Results of the self-consistent Hartree calculation for charge ordering with V0 = 30 meV.We show the net free energy gain for various charge patterns, including the stripe phase, and generalized Wigner crystal phases at n = 1/3, 1/4 and 1/7.On large 12 × 12 unit cells we also considered the possibility of amorphous charge states, which have the lowest free energy in some ranges of filling.

FIG. 3 .
FIG.3.DC resistivity in the vicinity of the n = 1/3 moiré-Wigner-Mott insulator (a-c) and for amorphous states close to n = 1/2 (d).Exactly at n = 1/3 (a) we find a strange metal phase at high temperatures, followed by activated behavior after charge-order sets in.For weak interactions, the Hubbard-Mott gap closes, and a Fermi liquid forms at low temperature.Heavy Fermi liquid behavior remains the dominant low-temperature behavior upon electron-doping (c); for hole-doping (a), Fermi liquid condensation occurs only at extremely low temperatures (not shown).In contrast, the amorphous state (d, shown for U/t = 38.3, at n = 1/2) is characterized by a weakly temperature-dependent resistivity turning from insulating to bad metallic behavior.The dashed line indicates the Mott-Ioffe-Regel limit.

FIG. 4 .
FIG.4.The density of states (DOS) as a function of energy for various doping and temperature.We only show here data for U/t = 19.1.At n = 1/3 (b), the spectrum is clearly split into a lower Hubbard band (LHB), a middle band consisting of the states on the unoccupied honeycomb sublattice, and the upper Hubbard band (UHB).The gap at the Fermi level is therefore a charge-transfer (CT) gap.Consequently, electron-doping (c) leads to the formation of a quasiparticle peak at the Fermi level whereas hole-doping (a) puts the Fermi level in the incoherent lower Hubbard band.Bad metal behavior, as shown in Fig.3, is due to a shifting chemical potential which can be clearly seen in the left figure.Upon increasing temperature the charge-transfer gap closes and a strange metal phase appears.The DOS for the amorphous state (d) has similar features (LHB, CT gap, UHB) but its features are smeared out due to the inhomogeneous charge patterns, leading to a soft gap in the spectrum and consequently bad metal behavior (Fig.3).

FIG. 5 .
FIG. 5. a.The Moiré potential in aligned WS2/WSe2 heterobilayers, based on the continuum model from Ref. [12] using our ab initio density functional calculations.There is a clear maximum at the W/Se2 stacking in the Moiré unit cell.b.The resulting bandstructure contains a flat band with bandwidth W = 15.8 meV separated from the other bands.c.The Wannier orbital corresponding to the flat band is approximately a Gaussian centered at the maximum of the Moiré potential.
APPENDIX D: PHASE DIAGRAM AND TRANSPORT PROPERTIES AROUND n = 1/3

FIG. 8 .
FIG. 8.The spatial variations in the LDOS over an area of 36 × 36 unit cells, for the states below (a-b), at (c) and above (d-f) the Fermi level.The spectroscopic mapping use the same color scale for all panels.T = 1.2K, n = 1/2, U/t = 38.3.

Fig. 10 FIG. 9 .
Fig.10shows the DC resistivity curves at n = 1/3 for various U across the transition.The MIT here has an essentially continuous character, similarly as for the Mott [1] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).

FIG. 10
FIG. 10.(a) DC resistivity curves for various U at n = 1/3; (b) the Arrhenius plot for the insulating curves in (a); dashed lins are obtained from linear fitting; (c) the activation gap extracted from (b); the red star ping points where metallic quasi-particles arise.