Excitonic Mott insulator in a Bose-Fermi-Hubbard system of moiré WS2/WSe2 heterobilayer

Understanding the Hubbard model is crucial for investigating various quantum many-body states and its fermionic and bosonic versions have been largely realized separately. Recently, transition metal dichalcogenides heterobilayers have emerged as a promising platform for simulating the rich physics of the Hubbard model. In this work, we explore the interplay between fermionic and bosonic populations, using a WS2/WSe2 heterobilayer device that hosts this hybrid particle density. We independently tune the fermionic and bosonic populations by electronic doping and optical injection of electron-hole pairs, respectively. This enables us to form strongly interacting excitons that are manifested in a large energy gap in the photoluminescence spectrum. The incompressibility of excitons is further corroborated by observing a suppression of exciton diffusion with increasing pump intensity, as opposed to the expected behavior of a weakly interacting gas of bosons, suggesting the formation of a bosonic Mott insulator. We explain our observations using a two-band model including phase space filling. Our system provides a controllable approach to the exploration of quantum many-body effects in the generalized Bose-Fermi-Hubbard model.


Introduction
The rich physics of the Hubbard model has brought fundamental insights to the study of manybody quantum physics [1].Initially proposed for electrons on a lattice, different fermionic and bosonic versions of this model have been simulated in various platforms, ranging from ultracold atoms [2] to superconducting circuits [3].Recently, bilayer transition metal dichalcogenides (TMDs) have become a versatile platform to study the Hubbard model thanks to the coexistence of several intriguing features such as the reduction of electron hopping due to the formation of moiré lattice with large lattice constant, and the existence of both intra-and inter-layer excitons.These characteristics have enabled the realization of numerous effects of many-body physics such as metal-to-Mott insulator transition [4][5][6][7][8][9], generalized Wigner crystals [10][11][12][13][14], exciton-polaritons with moiré-induced nonlinearity [15], stripe phases [16], light-induced ferromagnetism [17].Moreover, there have been recent exciting perspectives of exploring such effects in light-matter correlated systems [3,18,19].While typically the fermionic and bosonic versions of the Hubbard model are explored independently, combining these two models in a single system holds intriguing possibilities for studying mixed bosonic-fermionic correlated states [20,21].
In this work, we demonstrate Bose-Fermi-Hubbard physics in a TMD heterobilayer.We indepen- The population of the moiré lattice can be controlled via two independent parameters: the gate voltage changes the electronic filling factor (νe), and the optical pump creates a population of excitons, proportional to the input intensity.In the gray area, the system behaves as a mixed gas of bosonic and fermionic particles.As one approaches the upper limit (black line), the system becomes incompressible due to the saturation of single occupancy states.dently control the population of fermionic (electronic) particles by doping with a gate voltage (V g ), and the population of bosonic (excitonic) states by pumping with a pulsed optical drive of intensity I. Harnessing these two control methods, we realize strongly interacting excitons.In particular, we show the incompressibility of excitonic states near integer filling by observing an energy gap in photoluminescence, accompanied by an intensity saturation.Remarkably, we observe the suppression of diffusion, a strong indication of the formation of a bosonic Mott insulator of excitons.

Physical system and experimental design
To demonstrate these effects, we use a moiré lattice created by stacking two monolayers of WS 2 and WSe 2 , with symmetric top and bottom gates.Fig. 1a shows a schematic illustration of the heterobilayer device (See Supplementary Note 1 for details).Due to the type-II band alignment of the heterostructure (Fig. 1b), negative doping results in a population of electrons in the WS 2 subject to the moiré potential of the bilayer.The ratio between the density of this population and the density of moiré sites in the structure determines the so-called electronic filling factor (ν e ).The optical pump results in the formation of an energetically favorable interlayer exciton (X) [22], by pairing between an electron in WS 2 and a hole in WSe 2 (represented in Fig. 1b).In order to explore different regimes of Bose-Fermi-Hubbard model, we control the bosonic and fermionic populations by changing I and V g , respectively.This can be compared to ultra-cold atom implementation of Bose-Fermi mixture where the respective populations are fixed in each experiment [23].Before discussing our experimental observation, we discuss three limiting cases that determine the phase space of our system, as indicated in Fig. 1c.The corresponding physical scenarios are represented in panels d to f.First, in the weak excitation limit and low electronic filling factor (ν e ∼ 0) regime, the system's photoluminescence (PL) emission originates exclusively from the few X states in the quasi-empty lattice (panel d).This emission comes from excitons in lattice sites where they are the only occupant particles, namely, "single occupancy states" (X 1 ).Upon increasing ν e , the number of singly occupied sites decreases, and in the limiting case of ν e ≥ 1, as represented in panel e, the optically generated excitons can only form in lattice sites already occupied by charged particles.In this case, the required energy to form the exciton increases due to the on-site Coulomb repulsion, and hence the PL emission has new branch with higher energy than the previous regime.Consequently, the PL originates from lattice sites with an electron-exciton double occupancy (X 2 ).Finally, we consider the case where the electronic doping is below the threshold required to reach a fermionic Mott insulator (0 < ν e < 1) but I is strong enough to optically saturate the single occupancy states.The extra excitons create a number of sites with electron-exciton or exciton-exciton double occupancies (panel f).In this case, the PL emission corresponds to mixed contributions from exciton-exciton and exciton-electron interaction (U ex-ex and U ex-e ); the individual peaks cannot be distinguished in a single spectrum due to the broadness of linewidths.Therefore, in this regime, the emitted light is only a combination of the X 1 and X 2 PL emission.This interplay between exciton and electron occupancy can lead to situations in which the moiré lattice is completely filled with a mixed population of fermions and bosons, forming a hybrid incompressible state.Specifically, in the limit of weak electronic tunneling, excitons can form a Mott insulating state, in the remainder of sites that are not filled by electronic doping.Note the line in Fig. 1c denoting panel f is an asymptote since optical pumping can not fully saturate an exciton line.At ν e = 0, this intensity is denoted as I * .(See Supplementary Note 8 for details).

System's properties for varying electronic (fermionic) occupation
To experimentally investigate these regimes, we perform PL measurements, with varying pump power and backgate voltage.A detailed description of the optical setup can be found in Supplementary Note 2. We use pulsed excitation to achieve high exciton density while reducing thermal effects by keeping low average power.Experiments with CW excitation are consistent with the presented data, as shown in the Supplementary Note 6. Fig. 2a-c shows the PL dependence at three different intensities as schematically shown in panel d.Fig. 2a shows the normalized dopingdependent PL spectrum for low I (0.08µW/µm 2 ), which corresponds to low bosonic occupation.The fermionic occupation ν e is varied between 0 and 1.1.For low ν e , PL emission is detected only from X 1 .However, at V g ≈ 2.98V, we detect a transition in the PL emission to X 2 .This transition corresponds to the formation of X's in the presence of an incompressible fermionic Mott insulator [24,25].From the reflectivity measurement and calculations from a capacitor model, we attribute V g = 2.98 V to ν e = 1 (see Supplementary Note 3).The energy gap between X 1 and X 2 is ∆E ≈ 29meV, which corresponds to the on-site Coulomb repulsion energy between an electron and an exciton (U ex−e ).We elaborate on this energy gap later in the sub-section "Energy map along the phase space".The dim mid-gap features between X 1 and X 2 at ν e ∼ 1 are strongly position dependent and disappear at higher power.This indicates that such emission is from localized excitons.Fig. 2b shows the PL spectrum under pump intensity equal to 12.1 µW/µm 2 .It is worth noticing that the V g at which the PL signal from X 2 is detected, is lower than in panel a.The system is therefore in the regime depicted in Fig. 1f.Upon further increasing the pump intensity X 2 can be detected even at ν e = 0, as observed in Fig. 2c.In this case, the PL emission originates solely from double occupancy of excitons in a moiré lattice site, suggesting that, for high I, purely bosonic states of strongly interacting excitons are created.Comparing Fig. 2a and Fig. 2c, one can observe that in the former case, the emergence of the X 2 peak corresponds to a sharp suppression of X 1 , while in the latter case, both peaks coexist.This indicates the nature of the double occupancy: in the first scenario, the exciton is forming in the presence of an electron, and after its recombination, there are no other optical excitations in the system.In contrast, the coexistence of both peaks in panel c shows that upon double exciton occupancy, the recombination of X 2 precedes the recombination of X 1 .
From the observation described in the previous paragraph, we conclude that the detection of PL emission with X 1 and X 2 energies benchmarks the formation of exciton states in singly and doubly occupied lattice sites, respectively.At ν e = 0, the X 1 peak in Fig. 2c is blueshifted with respect to Fig. 2a.We associate this feature with a mean-field effect due to exciton-exciton interaction.As we increase the electronic doping, fewer sites are available to create X 1 excitons and on those occupied sites, only X 2 is created.Consequently, the effective population of X 1 excitons is decreased.Therefore the mean-field shift is suppressed to the point that at high filling (ν e ∼ 1) the X 1 energy is the same as in the case of low pump intensity.Next, in order to understand the interplay between fermionic and bosonic lattice occupancies in each regime, we perform a quantitative analysis of their respective integrated intensity.We extract these values from the collected PL spectra using a computational fitting method (see Supplementary Note 7 for further details).Fig. 2e-g displays this intensity dependence on ν e for the same I range of panels a-c.We notice that as electrons fill the system's phase space (upon increasing V g ), the number of accessible single-occupancy states decreases.As a consequence, the integrated intensity of X 1 reduces with increasing ν e .Remarkably, for each intensity, there is a critical ν e after which the PL emission of X 2 exceeds that of X 1 .The gate voltage at which the crossing takes place (V cr g ) is highlighted on each panel by a vertical dashed line.This line indicates a constant ratio between the X 1 and X 2 populations.The crossing takes place at lower ν e upon increasing I, as expected.In Fig. 2h, we track V cr g as a function of the total collected PL emission, which gives an indication of the total number of excitons in both X 1 and X 2 branches.We observe a clear trend: a higher total population of excitons results in a faster saturation of the single occupancy states and hence an increasing number of double occupancy states.

System's properties for varying excitonic (bosonic) occupation
Next, in order to trace the role of the optical pump and the optical saturation that leads to the formation of incompressible bosonic states, we investigate the PL for varying I for different ν e .In Fig. 3a-c, we focus on three different values of ν e , as indicated in panel d, and study the PL spectrum for increasing emitted PL power.For zero fermionic occupancy (panel a), X 2 contributes to the emission only at very high total PL emission intensity.In panels b and c, we increase the electronic doping to ν e = 0.7 and ν e = 0.95, respectively, and as expected, the total PL at which we detect X 2 decreases.In the low power region, panel c shows the PL emission from mid-gap states also observed in Fig. 2a.Apart from the energy gap in the emission, we observe a blueshift of the X 1 line with increasing PL power.Assuming the weak tunneling regime, this shift should be equal to U ex−ex ⟨x † x⟩, where x † is the creation operator of an exciton.For example, in Fig. 3b for total PL power at 2 counts/s, the bosonic occupation is ⟨x † x⟩ ≃ 0.2.This corroborates with the energy gap that occurs at 10 counts/s for an estimated unity filling (⟨x † x⟩ ≃ 1).We present a fully quantum theoretical analysis of this observation in Supplementary Note 10.Panels e-g show the intensities of X 1 and X 2 for the values of ν e in panels a-c.As expected, in panel e, one can observe that the intensity of the X 1 PL emission increases monotonically, and it starts to saturate only at very high total PL emission regimes.Upon filling the moiré lattice with one exciton or one electron per site, the X 1 PL intensity saturates.With higher ν e , the saturation occurs at lower I, as shown in panels f and g.Since this saturation corresponds to filling the single occupancy states, we associate it with the establishment of an incompressible bosonic Mott insulator.Note that this bosonic Mott insulator is in a drive-dissipative regime, similar to the demonstration in superconducting qubit systems [26].
To quantitatively analyze this saturation effect, we fit the X 1 PL power (P 1 ) to the function P 1 = P max 1 P P+Psat , where P is the total PL power.From the fitting, we extract P max 1 which is the asymptotic value of the X 1 emitted PL power, and P sat which determines the total PL of saturation.This functional form corresponds to the expected system behavior when the charge gap U is sufficiently large to permit the utilization of a phase-space filling model to treat both single and double occupancy states (details in Supplementary Note 8).Fig. 3f includes an example of the fitting function (dashed black line).According to our model, the value of P sat should decrease with increasing ν e because a lower excitonic population is required to achieve the incompressible states.The compiled data for the full range of ν e , shown in panel h with brown marks, is in good agreement with the expected trend.From the model, we can also infer that the quantity P sat /P max 1 should be independent of the electronic doping level because both quantities depend linearly on 1 − ν e ; higher electronic occupancy implies less single occupancy states available to host an exciton.The green marks in Fig. 3h represent this behavior, which is in good agreement with the model.We conclude that the saturation of single occupancy states is directly reflected in the intensity of X 1 , enabling the extraction of the conditions under which the incompressible states occur.Importantly, this enables a direct calibration of the bosonic and fermionic fractions in the system.

Exciton diffusion measurements
In order to further validate the incompressible nature of excitonic states, we perform diffusion measurements of the interlayer excitons [27].For a steady population of excitons created by a continuous-wave laser pump, the diffusion length carries information about the nature of the state: an incompressible bosonic state is expected to have a lower diffusion length than a weakly interacting one.We spatially image the diffusion pattern with spectral resolution and extract the diffusion length (L X1 ) of the single occupancy excitons.The choice of L X1 as an appropriate quantity to benchmark the incompressibility of bosonic Mott insulating states, assumes a constant exciton lifetime with varying population.This is supported by previous reports in the literature that show the independence of this quantity over three orders of magnitude of pumping power [28].The downward diffusion image has patterns that originate in the inhomogeneous surface of the bilayer.Although the inhomogeneities on that side hinder the extraction of L X1 , the optically-induced suppression of the diffusion length for constant ν e can be clearly observed in this region (Fig. 4a-b).
The population injected at y = 0 (dotted line) propagates, and its emission pattern is monitored along a range of 5 µm (dashed rectangle).The color scale is the same for both panels.Panel b shows a reduction of the diffused X 1 population in comparison to panel a.For the quantitative analysis of this observation, it is necessary to use a fitting routine, for which the smooth pattern on top of the injection point (y < 0) is more reliable.Figure 4c shows the extracted L X1 as a function of V g for different pump intensities from the exponentially decaying spatial diffusion pattern in this region.We provide more details about the analysis of the diffusion data in Supplementary Note 9.For low electronic density, the exciton diffusion length increases as the power is augmented.This trend, highlighted by the upward arrow, is in agreement with the expected behavior for weakly interacting bosons [28,29].Remarkably, as the electronic filling factor increases, the trend completely inverts (inset).This is a direct signature of the bosonic Mott insulator formation since the bulk is incompressible and the melting only occurs at the edge.

Energy map along the phase space
The implemented fitting algorithm allows us to track the changes in the energy of both species of excitons and the energy gap between them.These results are presented in Fig. 5. Panels a and b show the central energies of the peaks X 1 and X 2 in the space of parameters for which each peak is detectable.In the range where both of them can be detected, their energy difference ∆E (panel c) provides important information about the nature of the interactions taking place in the system.In the case of low electronic occupancy and high exciton density (top left corner of the panel), ∆E corresponds to the exciton-exciton interaction gap (Uex-ex ∼ 32 meV).Conversely, at high ν e and low exciton density (bottom right corner), this gap depends on the exciton-electron interaction (Uex-e ∼ 27 meV).The gradual change in the nature of the interactions taking place in the system along the parameters space is reflected in the change of ∆E.Interestingly, the largest energy gap takes place for states with high occupation of bosons and fermions (top right corner), which is consistent with a blueshift of the X 2 PL peak due to the high population of excitons with large Bohr radius repelling through dipolar interaction.

Discussion
In summary, we demonstrated a Mott insulating state of excitons in a hybrid Bose-Fermi Hubbard system formed in a TMD heterobilayer.While our incompressibility observation was based on spatially resolved diffusion in the steady-state limit, one can explore interesting non-equilibrium physics due to the relatively long lifetime of interlayer excitons.More generally, spatiotemporally resolved measurements, combined with independent tunability of fermionic and bosonic populations, make it possible to investigate the equilibrium and non-equilibrium physics of Bose-Fermi mixtures.Moreover, a quantum microscopic model capable of fully describing such a driven-dissipative Bose-Fermi mixture remains an open area of research.The novel experimental diffusion method used to benchmark the excitonic incompressibility opens exciting perspectives for the simulation of complex dynamics in many-body quantum systems that range from a single bosonic particle in a Fermi sea to a strongly interacting gas of bosons.Particularly intriguing examples are the optical investigation of charge and spin physics in integer and fractional fillings, e.g., Mott excitons [30? ] or spin liquids [31][32][33][34], and fractional Chern insulators [35,36].

Methods
Device fabrication: The WSe 2 /WS 2 heterostructure was fabricated using a dry-transfer method with a stamp made of a poly(bisphenol A carbonate) (PC) layer on polydimethylsiloxane (PDMS).All flakes were exfoliated from bulk crystals onto Si/SiO 2 (285 nm) and identified by their optical contrast.The top/bottom gates and TMD contact are made of few-layer graphene.The PC stamp and samples were heated to 60 • C during the pick-up steps and released from the stamp to the substrate at 180 • C. The PC residue on the device was removed in chloroform followed by a rinse in isopropyl alcohol and ozone clean.Sample transfer was performed in an argon-filled glovebox for improved interface quality.The electrodes consist of 3.5 nm of chromium and 70 nm of gold.They were fabricated using standard electron-beam lithography techniques and thermal evaporation.
Optical measurements: The sample is kept in a dilution refrigerator at a temperature of 3.5K.For PL measurements, we use a confocal microscopy setup.Our pumping source is a pulsed Ti:Sapphire laser tuned at 720nm (1.722eV), with a pulse duration of 100fs and a repetition rate of ∼ 80 MHz.Additionally, an optical chopper system at 800Hz is used to prevent sample heating while having a high pump intensity.The residual pump is removed with a spectral filter before collecting the PL emission in a spectrometer-CCD camera device.A complete description of the setup is presented in the Supplementary Note 2.
For the diffusion measurements, we used a continuous-wave (CW) laser.The rest of the optical measurement setup was similar.The WSe 2 /WS 2 heterostructure shown in Supplementary Figure S1 was fabricated using a dry-transfer method with a stamp made of a poly(bisphenol A carbonate) (PC) layer on polydimethylsiloxane (PDMS) [37].All flakes were exfoliated from bulk crystals onto Si/SiO 2 (285 nm) and identified by their optical contrast.The top/bottom gates and TMD contact are made of few-layer graphene.The PC stamp and samples were heated to 60 • C during the pick-up steps and released from the stamp to the substrate at 180 • C. The PC residue on the device was removed in chloroform followed by a rinse in isopropyl alcohol and ozone clean.Sample transfer was performed in an argonfilled glovebox for improved interface quality.The electrodes consist of 3.5 nm of chromium and 70 nm of gold.They were fabricated using standard electron-beam lithography techniques and thermal evaporation.

Supplementary Note 2. Optical measurements
The sample is kept in a dilution refrigerator at a temperature of 3.5K.For photoluminescence measurements, we use a confocal microscopy setup with an objective of magnification 70× and numerical aperture NA = 0.82.Our pumping source is a pulsed Ti:Sapphire laser tuned at 720nm (1.722eV), with a pulse duration of 100fs and a repetition rate of ∼ 80 MHz.A half-waveplate placed before a polarizing beam-splitter is rotated to control the pump power.Additionally, an optical chopper system at 800Hz is used to prevent sample heating while having a high pump intensity.A 750nm long-pass filter was used to block the residual pump laser before collecting the PL emission in a spectrometer equipped with a 300 grooves per mm diffraction grating and a CCD camera.The schematic of the setup is shown in Supplementary Figure S2.
For the diffusion measurements, we used a continuous-wave (CW) Ti:Sapphire laser tuned at 708nm.The rest of the optical measurement setup was similar.By applying a spatial filter to the reconstructed image of the diffusion pattern, we obtain spectral and spatial data of the exciton emission.This allows to image the spatial diffusion of each spectral component and

Supplementary Note 3. Electronic filling factor calibration
The electronic filling factor (ν e ) can be estimated by combining the information about the crystalline structure of the bilayer device and a parallel capacitor model, as follows: Determination of moiré density n 0 : In a purely fermionic Mott insulating state, the heterostructure will host one electron per moiré unit cell.The charge density in this case is given by n 0 , and it can be directly determined by the moiré periodicity through the relationship n 0 = 1/(L M 2 sin π/3).Here, L M = a/ √ δ 2 + θ 2 is the size of the moiré superlattice, δ = (a − a ′ )/a ≈ 4% is the lattice mismatch between WSe 2 (a = 0.328nm) and WS 2 (a ′ = 0.315nm), and θ is the twist angle between the two layers.Assuming 0 Determination of electron density n e : From a parallel capacitor model, we can deduce the expression for the electron density (n e ) in the bilayer.For a dual gate device, n e is given by: where V g is the symmetrically applied gate voltage and d t ≈ d b ≈ 40nm are the thicknesses of the top and bottom hBN dielectrics, respectively.They are determined from the optical contrast of the hBN flakes under the microscope.ϵ 0 is the vacuum permittivity, and ϵ r ≈ 3 is the relative permittivity of hBN [38].Having the moiré density n 0 and the electron density n e , the electronic (fermionic) filling factor can be expressed as ν e = n e /n 0 .From this model, we deduce that the gate voltage at which the electronic Mott insulator is established lies in the range 2.65 V< V g < 3.04 V.Where we are taking into account that the neutral region extends up to V g = 0.58 V.This estimation is in good agreement with our experimental data for reflectivity and PL.The results, displayed in Supplementary Figure S3 show a consistent change in the optical response of the device at V g = 2.98 V.For the reflectivity (panel a), we observe a shift in the intralayer exciton energy.
For the PL (panel b), a corresponding energy transition is observed for the interlayer exciton.Although this calibration holds for both the FIG.S3.Reflection and photoluminescence spectra of our WS2/WSe2 heterobilayer device.
hole-doped and the electron-doped sides, Supplementary Figure S3b shows an asymmetric behavior of the exciton interaction with the sign of the doping: the energy gap at ν e = 1 on the electron doping side is larger than the one in the hole doping side.This is in agreement with theoretical predictions [39] and experimental observations [20,40] that show a reduced value of the interaction strength due to the fact that holes and excitons reside in different high symmetry points of the moiré lattice.This observation indicates that electron doping is more suitable for the study of Fermi-Bose hybrid correlated states.

Supplementary Note 4. Determination of stacking angle via Second Harmonic Generation
We performed Second Harmonic Generation (SHG) measurements on our device in order to determine the stacking angle of the component monolayers.This is important because the system's properties vary for R-stacking and H-stacking.Using a 100 fs excitation at 850 nm with variable linear polarization, we are able to reconstruct the expected behavior of the SHG angular dependence for a crystalline structure with inversion symmetry breaking (Supplementary Figure S4a).We first identify spots on the sample with bilayer and monolayer regions by using the PL emission for the characterization (panel b).After identifying each spot, we pro-ceed to measure the SHG spectrum while keeping the fs laser power constant at ∼ 150µW .As demonstrated in ref. [41] and the Supplementary Material of ref. [16], the stacking angle can generate a constructive (destructive) interference in the SHG spectrum for 0°(60 • ) stacking.The strong suppression in the SHG efficiency observed in panel c, allows us to conclude that our sample corresponds to a bilayer with a 60 • stacking angle (AB or H stacking).

Supplementary Note 8. Modeling X1 and X2
To obtain further insights into the physical mechanism dominating the system dynamics, we make a theoretical analysis using a 2-band model, as mentioned in the main text.Our goal is to deduce an analytical expression for the observed saturation behavior of the PL power from the X 1 population.This is achieved by considering the steady state of the system in a continuous-wave excitation regime.A population of electron-hole pairs n eh is created by a Rabi-driving Ω(t).The charged plasma cascades down and relaxes to create n 1 singly occupied moiré sites (X 1 ) and n 2 doubly occupied sites (X 2 ).The formation of those states depends on the relaxation rate of the plasma (Γ relax ), the maximum number of available lattice sites (n max (ν e )), and a β factor, which determines the probability of the plasma to decay into the X 1 state.Under these conditions, the system's dynamics can be described by the set of equations: (S1c) where η is the efficiency of the optical generation of electron-hole pairs and Γ 1 (Γ 2 ) is the decay rate of the excitonic states in a singly (doubly) occupied site.In a steady state condition, from Eq. S1b, we can obtain an analytical expression for n 1 : In this regime n eh reaches an asymptotic value n eh = η Γ relax |Ω| 2 .To associate this expression with our experimental data, we take into account that n eh Γ relax is proportional to the total PL power and n 1 is proportional to the PL power (P 1 ) from the X 1 exciton band.After substituting variables we can conclude that the saturation behavior of X 1 should follow the functional form: where P max 1 = n max is proportional to the PL power from X 1 in an ideal Mott insulating state and P sat = Γ 1 n max /β determines the total PL emission at which X 1 saturates.This expression is used as a fitting model for the PL power from the X 1 states upon increasing total PL power, and provides a mathematical tool to estimate the bosonic occupation of the moiré lattice.From this theoretical treatment, one can notice that the bosonic Mott insulating states are reached asymptotically.This means that P 1 can be arbitrarily close to P max 1 upon a high vii enough total PL power.In other words, the line that determines the Mott insulating phases in Figs.1c, 2d, and 3d of the main text is mathematically not achievable.For this reason, the intensity denoted as I * in the referred figures, indicates the required pump intensity to drive the system into a state with a total population arbitrarily close to the lattice complete saturation.

Supplementary Note 9. Diffusion measurements
For the diffusion measurements, we use a continuous-wave laser tuned at 708 nm to create a steady population of excitons.By imaging the diffusion pattern with spectral resolution, we are capable of monitoring the diffusion properties of the X 1 and X 2 populations.We study the dependence of the diffusion length for each population as a function of V g and I.To improve the spatial resolution, the image of the diffusion pattern is magnified 200 times.We spectrally resolve the signals emitted from X 1 and X 2 .Next, we fit each spatial diffusion profile to a function of the form A exp (−x/L Xi )+b, where A is the PL power under the pumping laser spot, x is the propagation distance, L Xi is the diffusion length of X i and b is an offset that accounts for the base noise level.Supplementary Figure S8 (a-b) shows the obtained power dependence for L X1 at V g = 1.7 V and V g = 2.38 V. We observe that L X1 increases with pump power at 1.7 V but decreases at 2.38 V.This inversion in the power dependence accounts for the formation of incompressible excitonic states.It is worth mentioning that the scale of the I axis is not linear, because the power was modified by changing the polarization of the pump laser with a half-wave plate.Panel c shows the fitting subroutine used to extract the diffusion length of X 1 .Panel d shows the measured value of L X1 for a range of V g to highlight the inversion of the intensity dependence.The same analysis was performed for the diffusion of the X 2 states.The data (Supplementary Figure S9a) shows a shorter diffusion length, which makes it cumbersome to detect changes in L X2 due to the diffraction limit.However, in contrast to L X1 , the L X2 increases with increasing V g , except for the reduction at ν e ∼ 1, which can be attributed to the behavior of a population of excitons in presence of a fermionic Mott insulating state.
We finally analyze the diffusion data with no spectral resolution.For ν e < 1, one can observe that the suppression of the diffusion takes place for lower V g upon increasing pump intensity; a consequence of the bosonic saturation of the lattice due to the optical pump.
Delving deeper into the spatial diffusion, we identify a transfer of population from X 2 to X 1 as the excitons move away from the pump region.As observed in Supplementary Figure S10, for high pump power, the population of X 1 away from the excitation spot is larger than under it, a manifestation of the mentioned effect.At the excitation spot, the high exciton density leads to the formation of a localized Mott insulating state and therefore a large population of X 2 .However, as the distance from the pumping region increases, part of the population transfers from X 2 to X 1 .Notice how, for high excitation intensity, the X 1 signal becomes stronger away from the injection spot than under it.In this case, the X 2 population is acting as a reservoir for X 1 along the diffusion path.Importantly, this behavior does not change the interpretation of our results, but constitutes proof of the robustness of the effect; the population transfer from X 2 to X 1 would favor an apparently longer propagation constant for X 1 , which, at most, would hinder the measured suppression of the diffusion.creasing excitonic population (Supplementary Figure 3a of main text): a blueshift at low occupancy and the emergence of a PL peak from a gapped state after a threshold PL power.In this "reduced" system, where the population is fully bosonic, we can model the system and account for the blueshift and the spectral gap by using the following Lindbladian master equation: where Ĥ corresponds to the Bose-Hubbard Hamiltonian: and the operators x(x † ) stand for the annihilation (creation) of an exciton particle, ω x is their energy, and U ex−ex is the on-site particle repulsion interaction energy.We consider two Lindbladian terms to take into account the laser pump and exciton loss.Specifically, the operator for channel n is: that includes the jump operators accounting for two incoherent processes: laser pump and exciton losses, with operators Ĉp and Ĉl , respectively: Using the quantum regression theorem for the expected value ⟨x † x⟩, we calculate the spectral The results, displayed in Supplementary Figure S11 show a very good agreement with the experimental results.This figure provides important insights into the nature of U ex−ex .It shows that the on-site exciton repulsion leads to a spectral blueshift even before the system reaches an insulating state.Importantly, such a master equation cannot explain the intensitydependent behavior of X 1 .A unified model including higher-order jump operators (and hence more fitting parameters) may be required.We anticipate this master equation description and the rate equations presented in section VIII being simplifications of a unified theory, which needs further investigation.tively.Supplementary Figure S16a presents the X 1 and X 2 energy evolution of device D2 under a non-resonant pump (633 nm), where we can observe a strong photo-doping effect -a reduction of X 1 intensity as a function of pumping power, consistent to the gate voltage-dependent PL spectrum at low power (Fig. 2e in the main text).However, while using the resonant pump (720 nm) for D1, as indicated in panel b, the X 1 emission monotonically increases for the whole range of pump power, an observation confirms that no photo-doping effect exists in the system.

FIG. 1 .
FIG. 1.(a) Schematic of the WS2/WSe2 dual-gate device.The TMD heterobilayer is embedded between two symmetric gates: top gate (TG) and bottom gate (BG).(b) Depiction of the type-II band alignment of the bilayer.The blue and red curves denote bands from WS2 and WSe2, respectively.The shaded ellipse indicates the formation of interlayer excitons composed of an electron from the WS2 conduction band and a hole from the WSe2 valence band.(c) Phase diagram of the system.The population of the moiré lattice can be controlled via two independent parameters: the gate voltage changes the electronic filling factor (νe), and the optical pump creates a population of excitons, proportional to the input intensity.In the gray area, the system behaves as a mixed gas of bosonic and fermionic particles.As one approaches the upper limit (black line), the system becomes incompressible due to the saturation of single occupancy states.(d-f) Interlayer exciton formation under optical excitation for three different regimes governed by the pump intensity (I) and νe: (c) low I and νe ∼ 0, (d) low I and νe ∼ 1, (e) high I and 0 < νe < 1. X1 (X2) denotes PL emission from singly (doubly) occupied moiré lattice sites.X2 can originate from either electron-exciton (Uex−e) or exciton-exciton (Uex−ex) double occupancies.
FIG. 1.(a) Schematic of the WS2/WSe2 dual-gate device.The TMD heterobilayer is embedded between two symmetric gates: top gate (TG) and bottom gate (BG).(b) Depiction of the type-II band alignment of the bilayer.The blue and red curves denote bands from WS2 and WSe2, respectively.The shaded ellipse indicates the formation of interlayer excitons composed of an electron from the WS2 conduction band and a hole from the WSe2 valence band.(c) Phase diagram of the system.The population of the moiré lattice can be controlled via two independent parameters: the gate voltage changes the electronic filling factor (νe), and the optical pump creates a population of excitons, proportional to the input intensity.In the gray area, the system behaves as a mixed gas of bosonic and fermionic particles.As one approaches the upper limit (black line), the system becomes incompressible due to the saturation of single occupancy states.(d-f) Interlayer exciton formation under optical excitation for three different regimes governed by the pump intensity (I) and νe: (c) low I and νe ∼ 0, (d) low I and νe ∼ 1, (e) high I and 0 < νe < 1. X1 (X2) denotes PL emission from singly (doubly) occupied moiré lattice sites.X2 can originate from either electron-exciton (Uex−e) or exciton-exciton (Uex−ex) double occupancies.

FIG. 2 .
FIG. 2. (a-c) Normalized PL spectrum as a function of gate voltage (νe) for three different pump intensities: I = 0.08 µWµm 2 (a), I = 12.1µW/µm 2 (b) and I = 1229µW/µm 2 (c).The peaks associated with single (X1) and double (X2) occupancy are indicated on each panel.The dashed lines indicate the gate voltages at which the PL intensity X2 exceeds X1.The dashed black lines of panel (d) indicate the measurement ranges of (a-c).(e-f) Evolution of the PL intensity for X1 (red) and X2 (blue) as a function of gate voltage for the same values of pump intensities displayed in panels (a-c).The electron filling factor at which X2 exceeds X1 decreases as pump intensity increases.Panel (h) shows the gate voltage at which the intensity of X2 exceeds that of X1, as a function of the total PL intensity.The error bars represent the standard errors for the parameter estimates in the fitting routine.

FIG. 3 . 1 (
FIG. 3. (a-c) Normalized PL spectrum as a function of the total collected PL power for three different electronic filling factors.The peaks associated with single (X1) and double (X2) occupancy are indicated on each panel.Panel (d) indicates the ranges of I and νe for the measurements shown in panels (a-c).(e-f) evolution of the PL power for X1 (red) and X2 (blue) as a function of the total collected PL power for the same values of νe displayed in panels (a) to (c).Panel (f) displays the fitting function (dashed black line) employed to extract Psat and P max 1 (described in the text).(h) Evolution of Psat (brown) as a function of the gate voltage (νe).As expected from our phase-space filling model, its value reduces with increasing filling factor.The quantity Psat/P max 1 (green) shows good agreement with the theoretical model.The error bars represent the standard errors for the parameter estimates in the fitting routine.

FIG. 4 .
FIG. 4. Spectrally and spatially-resolved diffusion pattern at νe = 0.73 (Vg = 2.34V) for low (a) and high (b) I.The dashed rectangle highlights the region where the suppression of diffusion can be observed.(c) Exciton diffusion length as a function of the gate voltage for a range of νe and for different input intensities.For low νe the diffusion length increases with I due to exciton repulsive interaction.Upon further filling the moiré lattice, the trend inverts, indicating the optical realization of incompressible states.The inset is a zoom-in of the red dotted rectangle to highlight the reduction of LX1 with increasing I.The error bars represent the standard errors for the diffusion length estimated from the exponential fitting.

FIG. 5 .
FIG. 5. Energy of the X1 (a) and X2 (b) PL emission as a function of gate voltage and pump intensities.The white areas correspond to the range of parameters where the corresponding peak completely vanishes.When X1 and X2 coexist, we extract the energy difference, as shown in panel (c).

3 .
Contributions B. G., D.G.S.F, S.S., and M.H. conceived and designed the experiments.K.W. and T.T. supplied necessary material for the fabrication of the sample.B.G., D.S., and R.N. designed and fabricated the sample.J.V. and S.M. collaborated with the preparation of the setup at its initial stage.B.G., D.G.S.F., and S.S. performed the experiments.B.G., D.G.S.F., S.S., T.S.H., M.J.M., M.X., A.I., Y.Z., and M.H. analyzed the data and interpreted the results.T.S.H. and M.H. elaborated on the theoretical models presented in the manuscript.B.G., D.G.S.F., S.S., M.J.M., and M.H. wrote the manuscript, with input from all authors.Electronic filling factor calibration 4. Determination of stacking angle via Second Harmonic Generation 5. Total emitted PL vs. Pump Intensity 6. Pulsed vs. CW excitation 7. Data analysis 8. Modeling X 1 and X 2 9. Diffusion measurements 10.System's spectrum in the purely bosonic limit 11.Reproduction of the observations in Device D2 12.Comparison of U ex−ex and U ex−e with previously reported values 13.Optical signatures of photo-doping effect Supplementary Note 1. Device fabrication FIG. S1.Microscope image of the WS2/WSe2 heterostructure device.
FIG. S2.Schematic of the experimental setup.HWP=Half-wave plate, BS=Beam splitter FIG. S7.A typical PL spectrum (Vg = 2.1 V and I = 2226 µW/µm 2 ) with the fitting functions obtained from the numerical method.Inset shows a wide PL spectrum at low and high pump powers with no WSe2 intralayer excitons or trions.
FIG. S8. (a-b) Pump intensity dependence of the X1 spatial diffusion for νe ∼ 0.47 (Vg = 1.7V ) (a) and νe ∼ 0.75 (Vg = 2.38V ) (b).The panels show the inversion of the intensity dependence: from a diluted gas of bosons in panel a to a bosonic Mott insulator in panel b.(c) Exponential decay fitting of a typical X1 diffusion pattern.(d) LX1 for a reduced range of Vg to highlight the power dependence inversion.The error bars represent the standard errors for the diffusion length estimated from the exponential fitting.
FIG. S9.(a) LX1 as a function of Vg for different pump intensities.In this panel, we include dashed lines corresponding to the crossing gate voltage at each power.They are a guide for the eye of where the incompressibility is expected to manifest at each power.(b)LX2 as a function of the gate voltage for a range of νe and for different I (c) LX as a function of the gate voltage for a range of νe and different I.The error bars represent the standard errors for the diffusion length estimated from the exponential fitting.

∂
FIG. S10.Diffusion of the system for νe = 0 and two different pumping intensities.a) 9.42 µW/µm 2 and c) 1270 µW/µm 2 .The right panels (b and d) show the horizontally integrated intensity.The transfer of population from X2 to X1 becomes evident in panel d, where the X1 intensity is larger outside the pumping spot due to a large reservoir of X2 excitons that decay into X1 as they propagate.
FIG. S11.Continuous line: Main peak position of the calculated spectral function for a Bose-Hubbard model under the Master equation formalism.The blue circles correspond to the obtained experimental data (from Fig. 3a of the main text) for the PL power dependence of the spectrum at νe = 0. x FIG. S12.(a) Gate voltage-dependent normalized PL in device D2 for a pump intensity of 0.048 µW/µm 2 .The splitting between X1 and X2 (Uex−e) is 34 meV.Inset shows an optical image of device D2.(b) Power-dependent normalized PL in device D2 at νe = 0.The splitting between X1 and X2 (Uex−ex) is 35 meV.