Nonlinear spectroscopy of semiconductor moiré materials

Abstract

We use time-resolved nonlinear pump–probe measurements to reveal features of semiconductor moiré materials not accessible to linear spectroscopy. With an intense, red-detuned pump pulse, we generate a high density of virtual excitons or exciton–polarons in various moiré minibands. A broadband probe pulse in turn measures the response of all optical resonances induced by the pump-generated excitations. We generically observe a coherent blue shift originating from contact-like exciton–exciton interactions. At charge neutrality, these measurements allow us to assess the spatial overlap between different optical excitations and to observe signatures of a bound biexciton state between two different moiré exciton modes. In contrast to electron doped monolayers, spatially confined moiré attractive polarons behave as an ensemble of non-interacting two-level emitters, exhibiting an electron-density-independent ac-Stark effect. Tuning the pump laser into resonance with the attractive polaron, we demonstrate the filling of the moiré lattice with localized polarons and thereby realize a nonequilibrium Bose–Fermi mixture in moiré flat bands.

I Introduction

Semiconductor moiré materials have emerged as a rich playground for exploration of strongly correlated electrons [1, 2, 3]. Recent experiments based on twisted bilayers of transition metal dichalcogenides (TMDs) have enabled the observation of a wealth of many-body states, ranging from correlated Mott–Wigner states [4, 5, 6], through kinetic magnetism [7] to fractional Chern insulators [8, 9]. These discoveries were made using various experimental techniques, such as transport, STM spectroscopy, or linear optical spectroscopy [10]. The link between the optical measurements and the electronic state stems from the dynamical dressing of tightly bound TMD excitons by itinerant charges, leading to the formation of attractive and repulsive exciton–polarons  [11, 12, 13, 14, 15, 16]. Since exciton–polarons are sensitive to both the charge and magnetic order of the electron system, they provide unequivocal signatures of correlated electrons  [17, 18, 7]. A common feature of all the measurement techniques implemented to date is that they exclusively measured linear response of electrons to external fields.

In this work, we use nonlinear optical pump–probe measurements to investigate features of semiconductor moiré materials not accessible through linear spectroscopy. In particular, we investigate their optical excitations, which emerge from a complex interplay between the moiré potential and the interaction bewteen excitons and a many-body electronic state in a flat band. Our approach relies on an intense red-detuned pump beam to generate a significant density of virtual excitations, and on a weak broadband probe pulse to monitor the subsequent modification of the optical spectrum. By detuning the pump away from any resonances, we alleviate real-absorption-induced modification of the electronic state and look at the coherent scattering response of the system. When the heterobilayer is charge neutral, we observe a blue shift of moiré exciton resonances in co-circularly-polarized pump-probe measurements [19, 20, 21, 22, 23, 24, 25, 26]: these measurements reveal the extent of spatial overlap between distinct exciton modes. We find further confirmation of this overlap under cross-circularly-polarized excitations [27, 28, 29, 30, 31] where we identify biexciton Feshbach resonances associated with bound states of both the same and different moiré exciton modes [32]. To further study the nature of these biexciton states, we develop a theoretical model for the scattering of excitons in a moiré and find good qualitative agreement with the experiment. When the moiré material is electron doped, nonlinear pump-probe measurements reveal an electron-density-independent blue shift of the attractive polaron resonance that scales linearly with the inverse pump-laser-detuning: such a single-emitter-like ac-Stark shift [33, 34, 35] demonstrates that the attractive polaron resonance should be considered as an ensemble of spatially localized and non-overlapping independent trion excitations.

Figure 1: a Schematic of the device. The TMDs MoSe₂ and WS₂ form the moiré material that we investigate. It is encapsulated in two $\approx 35\,$nm thick hBN flakes, and dual gated with graphite electrodes to enable an independent control of the chemical potential and electric field. The latter is kept close to zero, such that only MoSe₂ is doped at low densities, given the band alignment (b). The evolution of the reflection spectrum of MoSe₂ as a function of the electronic density (c) reveals several moiré exciton and polaron resonances, which we investigate. We rely on a pump–probe scheme (d) where an intense red-detuned laser generates a virtual population of moiré excitations. The interaction between this background and a test excitation generated with a probe laser are then measured.

Our measurements are carried out in a $\simeq 0^{\circ}$-degree stacked MoSe₂/WS₂ heterostructure, exhibiting a Type I band alignment where the lowest (highest) energy moiré conduction (valence) band resides in MoSe₂ (Fig. 1 a, b). The electron density dependent reflection contrast exhibits four bright resonances, which we identify as MX₁, MX₂, AP and MX₃ to be consistent with the notation used in an earlier publication [7] (Fig. 1 c). A strong pump laser with a finite red detuning from a given resonance generates a large virtual population of the corresponding moiré exciton species/modes that exists only during the $\approx 0.2$ ps duration of the pump pulse. We estimate a maximum virtual exciton population in MX₁ of $\approx 10^{12}$ cm^-2 for a pump detuning of $\delta_{1}\approx 20$ meV (for more details, see the Appendix). A weak, broadband probe pulse then measures the energy shift of all excitonic species concurrently (Fig. 1 d). Throughout this work, we set the pump laser detuning from the excitonic resonances to be much smaller than the exciton binding energy: In this limit, the dominant contribution to the light shift of the resonances originates from exciton–exciton interactions [19, 20, 21, 22, 23, 24, 25]. We determine both intra- and inter-species interaction strengths between the moiré exciton or polaron modes as function of the electron density by measuring the light shift in this small detuning limit.

II Interactions between moiré excitons

We begin by performing non-linear spectroscopy at charge neutrality in order to explore exciton–exciton interactions in the presence of a moiré. When the material is free of itinerant electrons, two bright moiré exciton resonances MX₁ and MX₂ are visible in the normalized reflection spectrum of MoSe₂, $\Delta R=(R-R_{\mathrm{bg}})/R_{\mathrm{bg}}$ with $R$ and $R_{\mathrm{bg}}$ the reflection spectrum from the sample and background, respectively (see Fig. 1c). To explore the interactions involving these two resonances we consider two different scenarios: co- and cross-circular polarized light of the pump and probe beams. These scenarios allow us to either probe interactions between excitons in the same valley (co-circular) or in opposite valleys (cross-circular).

II.1 Co-circularly polarized light

Figure 2: Light shift of moiré excitons for co-circularly polarized pump and probe lasers. (a) Schematic of the moiré exciton energy level at charge neutrality. (b) Reflection spectrum as a function of the pump–probe delay, showing a clear blue shift of the lowest and brightest resonance MX₁ at zero time delay. The differential spectrum (c), obtained by subtracting from (b) a reference spectrum (acquired at $\tau<-2.5\,$ps), reveals a blue shift of the upper and darker resonance MX₂ as well. (d) The dependence of the blue shift of MX₁ (blue circles) on the pump detuning scales as $1/\delta_{1}^{2}\propto n_{1}$ (blue line), and stems from MX₁–MX₁ interactions. The shift of MX₂ (red squares), is well fitted either by a $\propto 1/\delta_{1}^{2}$ dependence (dashed red line) or by $\sqrt{n_{1}n_{2}}\propto 1/(\delta_{1}\delta_{2})$ (dotted red line). While we are not able to precisely deconvolve the role of these two potential contributions, both imply significant MX₁-MX₂ interaction and hence spatial overlap between these two bright moiré excitons.

Beginning with the scenario of co-circular polarization, we consider the evolution of $\Delta R$ as a function of the delay $\tau=t_{\rm probe}-t_{\rm pump}$, between the pump and probe pulse, as shown in Fig. 2 b (here the pump pulse is red-detuned from the MX₁ resonance by $\delta_{1}=20$ meV). The dominant coupling mechanism between two same-valley excitons are electron- and hole-exchange interactions [36, 37], which leads to short-range repulsion. Correspondingly, we observe a clear blue shift of the brightest exciton MX₁ for $|\tau|\lessapprox 0.2\,$ps, that is, when the two pulses overlap in time (for more detail on the experimental setup see the Appendix). To better assess the pump-induced modifications of weaker resonances, we use the differential reflection spectrum $\Delta R(t)-\Delta R_{\rm ref}$, where $\Delta R_{\rm ref}$ is a reference spectrum obtained when the probe pulse hits the sample significantly before the pump ($\tau\lessapprox-2.5\,$ps): Figure 2 c shows that a smaller blue shift of MX₂ is discernible in $\Delta R(t)-\Delta R_{\rm ref}$. In a mean-field picture (detailed in App. H.6), the light shift $\Delta_{i}$ of MX_i ($i=1,2$) can be expressed as

$$\Delta_{1}=2u_{11}n_{1}+2u_{12}n_{2}+4k_{1}\sqrt{n_{1}n_{2}}\,,$$ (1a)

$$\Delta_{2}=2u_{22}n_{2}+2u_{12}n_{1}+4k_{2}\sqrt{n_{1}n_{2}}\,,$$ (1b)

where $n_{i}$ is the density of MX_i excitons and we have introduced four different interaction terms: $u_{i,j}$ correspond to MX_i+MX${}_{j}\leftrightharpoons$ MX_i+MX_j, while $k_{i}$ correspond to MX₁+MX${}_{2}\leftrightharpoons$ MX_i+MX_i. The density scales as $n_{i}\propto 1/\delta_{i}^{2}$ with the detuning $\delta_{i}$ of the pump from the MX_i resonance. Therefore, by tuning the pump laser frequency, one can change the density imbalance between the two excitons, and in principle deconvolve the contribution of each scattering processes to the light shift. In practice, we observe signatures of an incoherent response for a blue-detuned pump, so we focus exclusively on red detunings. We consequently always have $\delta_{1}<\delta_{2}$ and thus $n_{1}>n_{2}$ (see Fig. 2 a). This imbalance is further amplified by the oscillator strength difference between MX₁ and MX₂.

Figure 2 d shows the light shifts $\Delta_{1,2}$ as a function of $\delta_{1}$ in a range where $n_{1}\gtrsim 10\,n_{2}$. As a result, the light shift of MX₁ is dominated by MX₁–MX₁ interactions. From a fit we obtain the interaction strength $u_{1,1}$, which we find to be larger by a factor $u_{1,1}/u_{\rm ex}\approx 1.6\pm 0.2$ than the moiré-free exciton–exciton interaction strength $u_{\rm ex}$, measured on a monolayer MoSe₂ region of the same device. We tentatively attribute this enhancement to the reduced spatial extent of the MX₁ center-of-mass wavefunction within the moiré unit cell, which in turn increases the overlap between MX₁ excitons for a given average density.

Remarkably and despite of the large imbalance $n_{2}\ll n_{1}$, we observe a substantial light shift of MX₂, strongly increasing as the pump wavelength approaches the MX₁ resonance. This shift $\Delta_{2}$ is not well reproduced by a $1/\delta_{2}^{2}$ dependence, and attests to substantial inter-species interactions. A fit to our data does not enable us to precisely disentangle the respective contributions of the MX₂+MX${}_{1}\leftrightharpoons$ MX₂+MX₁ and MX₂+MX${}_{1}\leftrightharpoons$ 2MX₁ processes. Nevertheless, we estimate that $0.3\lesssim u_{1,2}/u_{1,1}\lesssim 0.8$ while we are not able to reliably estimate the $k_{i}$ parameters (see App. C for more details on the fitting procedure).

These observations suggest a significant spatial overlap of MX₁ and MX₂ excitons and invalidates a simplistic picture of moiré exciton modes that are tightly confined around different high-symmetry points of the moiré potential. Moreover, our experiments are in reasonable qualitative agreement with the findings of moiré exciton wavefunctions obtained for the same structure using the continuum model [38]. In particular, we find that the continuum model predicts $u_{1,1}/u_{\rm ex}=2.2$ and the ratio of interactions to be $u_{1,2}/u_{1,1}=0.6$ (see Appendix H for details).

We point out that Eq. (1) can be understood as arising from first order perturbation theory. This approach is valid in the limit for small enough pump power/large enough pump detuning. We find an empirical confirmation that we are indeed working in that regime by observing a linear dependence of the light shifts $\Delta_{i}$ with the pump intensity (App. B). Nevertheless, we can envision interesting higher order effects arising from the mixing of the excitonic states. In particular, mixing of optically bright and dark excitons, could be observed as the emergence of new resonances in the reflection spectrum.

Figure 3: Coupling to biexcitonic states. (a) Reflection spectrum as a function of the pump–probe delay when the pump is red-detuned by $22\,$meV from MX₁ and the pump and probe laser are cross-circularly polarized. The differential spectrum (b), obtained by subtracting from (a) a reference spectrum (acquired at $t<-2.5\,$ps), reveals a blue shift of both resonances MX₁ and MX₂. A similar measurement carried out for a detuning of $45\,$meV from MX₁ shows instead a red shift of both resonances (c, d). The dependence of the light shifts of MX₁ (blue dots) and MX₂ (red dots) displays a sign change for a pump detuning $\delta_{1}\approx 36\,$meV and $\delta_{1}\approx 28\,$meV, respectively (e). (f) The corresponding theory simulation for the light-shifts of MX₁ (blue) and MX₂ (red) based on DFT moiré parameters. (g) The calculated spectrum which is shown in three different energy sectors with zero (defined as zero energy), one and two excitons. In the one exciton sector there are the optically bright MX₁ and MX₂ excitons. In the two exciton sector (focusing on zero centre-of-mass quasi-momentum) there exists a band of MX₁-MX₁ excitons (light blue) as well as other bands (gray), and bound biexcitons (green lines). We find that the sign change in the light shift of MX₁ and MX₂ in (f) can be attributed to the blue and red arrows (indicating exciton–biexciton transitions), respectively.

II.2 Cross-circularly polarized light

Since electron and hole exchange interactions are suppressed for excitons generated in opposite valleys, one may naively assume that the cross-circularly polarized scenario should yield a significantly smaller light shift. However, the bare interaction of opposite valley excitons is attractive, and supports a bound biexciton state, which has significant implications.

In order to understand how the biexciton state affects these interactions, we briefly review the simpler scenario in which there is no moiré (such as in the case of monolayer TMDs). In this case it is known that the presence of the biexciton resonance leads to an additional contribution to the light shift, scaling linearly with the inverse of the two-photon detuning $\delta_{\rm b}^{-1}=(-E_{\rm b}+\delta_{\rm ex})^{-1}$, where $\delta_{\rm ex}$ is the detuning of the pump from the exciton resonance, $E_{b}$ is the biexciton binding energy and the probe is assumed resonant with the exciton. A hallmark of the biexciton is the ac Stark effect, where the light shift changes in sign for detunings in the vicinity of $E_{\rm b}$ [27, 28, 29, 30, 31]. The two-photon resonance condition $\delta_{\rm b}=0$ can be considered as a biexciton Feshbach resonance where the effective interactions between the pump and probe generated excitons changes from being attractive ($\delta_{\rm b}>0$) to repulsive ($\delta_{\rm b}<0$).

We now extend these concepts to the case of excitons in moiré materials. By measuring the ac Stark effect of MX₁ and MX₂ excitons under cross-circularly-polarized pump-probe lasers as a function of $\delta_{1}$, we explore the nature of biexciton resonances in moiré materials: while a monolayer hosts a single zero momentum biexciton, we find that the presence of moiré potential leads to multiple zero quasi-momentum biexcitons. We first observe that both the MX₁ and MX₂ excitons blue shift close to resonance ($\delta_{1}\approx 22$ meV, Fig. 3 a and b) and red shift at larger detunings ($\delta_{1}\approx 45$ meV, Fig. 3 c and d). Figure 3 e shows the light shift as a function of the pump detuning: from a heuristic fit with the function $\Delta_{i}=a\tan^{-1}{\left(\frac{\delta_{1}-E_{\mathrm{b},i}}{b}\right)}+c$ (with $i=1,2$, $a,b,c$ fitting parameters), we extract the energy of two biexciton states $E_{\mathrm{b,1}}$ and $E_{\mathrm{b,2}}$ determining the sign change in the ac Stark shifts of MX₁ and MX₂. From the fit to the MX₁ light shift, we extract $E_{\rm b,1}\approx 36\,$meV, which suggests that the ground state moiré biexciton is more strongly bound than the previously measured MoSe₂ monolayer biexciton (with $E_{\rm b}\approx 28\,$meV [25]). The fit to the MX₂ light shift shows the presence of a second biexciton state, with a mixed MX₁/MX₂ character, and a binding energy $E_{\rm b,2}\approx 28\,$meV measured with respect to an unbound MX₁ and MX₂ excitons. Importantly, $E_{\rm b,2}$ is smaller than the energy splitting between MX₁ and MX₂ ($\approx 40$ meV), and this second biexciton is higher in energy than two unbound MX₁ excitons. This should be contrasted with the monolayer scenario, where the bound state by definition has lower energy than two unbound zero-momentum excitons.

II.3 Theoretical model

To further study the nature of the biexciton states, we employ a theoretical model for the scattering of excitons in a moiré potential, which draws inspiration from recent treatments of atoms interacting in low-dimensional optical lattices [39]. The Hamiltonian we consider is

$$\hat{H}=\sum_{\mathbf{k},\sigma,\lambda}E_{\mathbf{k}\lambda}\hat{X}^{\dagger}_{\mathbf{k}\lambda\sigma}\hat{X}_{\mathbf{k}\lambda\sigma}$$

$$+\sum_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}\\ \lambda_{\uparrow},\lambda_{\uparrow}^{\prime},\lambda_{\downarrow},\lambda_{\downarrow}^{\prime}\end{subarray}}\mathcal{V}^{\lambda_{\uparrow}^{\prime}\lambda_{\downarrow}^{\prime}}_{\lambda_{\uparrow}\lambda_{\downarrow}}(\mathbf{k},\mathbf{k}^{\prime},\mathbf{q})\hat{X}^{\dagger}_{\mathbf{k}-\mathbf{q}\lambda_{\uparrow}^{\prime}\uparrow}\hat{X}^{\dagger}_{\mathbf{k}^{\prime}+\mathbf{q},\lambda_{\downarrow}^{\prime}\downarrow}\hat{X}_{\mathbf{k}^{\prime}\lambda_{\downarrow}\downarrow}\hat{X}_{\mathbf{k}\lambda_{\uparrow}\uparrow},$$ (2)

where $\hat{X}_{\mathbf{k}\lambda\sigma}$ annihilates a moiré exciton with energy $E_{\mathbf{k}\lambda}$, quasi-momentum $\mathbf{k}$ and band index $\lambda$. Here the valley of the exciton is indexed by the spin of the electron inside the exciton, i.e., $\sigma=\uparrow,\downarrow$. The interaction term is written generically, and represents a short-range interaction potential in the basis of the moiré excitonic states: the potential is fixed by enforcing that it yields the correct biexciton binding in the absence of moiré. While this has been measured in monolayer as $E_{\rm b}\approx 28$ meV [25], we reduce it to 20 meV since we anticipate screening from the WSe₂ layer. The parameters for the moiré potential — which uses the continuum approximation — are derived from large-scale Density Functional Theory (DFT) [7]. We find that in between MX₁ ($\lambda=0$) and MX₂ resonances ($\lambda=3$), there are two dark excitonic resonances ($\lambda=1;2$). While we relegate most details of our model — including a detailed description of the interaction potential and the moiré parameters — to Appendix H, we now provide an overview of the model, including its comparison to experiment.

The scattering physics between an $\uparrow$ and $\downarrow$ exciton is fully characterised by the $T$ matrix, which corresponds to an infinite sum over exciton-exciton scattering events. We can study the experimentally relevant zero quasi-momentum biexcitonic states by calculating the matrix element of the $T$ matrix given by

$$T^{\lambda_{\uparrow}\lambda_{\downarrow};\mathbf{0}}_{\lambda_{\uparrow}^{\prime}\lambda^{\prime}_{\downarrow};\mathbf{0}}(E)\equiv\bra{0}\hat{X}_{\mathbf{0}\lambda_{\downarrow}\downarrow}\hat{X}_{\mathbf{0}\lambda_{\uparrow}\uparrow}\hat{T}(E)\hat{X}^{\dagger}_{\mathbf{0}\lambda_{\uparrow}^{\prime}\uparrow}\hat{X}^{\dagger}_{\mathbf{0}\lambda_{\downarrow}^{\prime}\downarrow}\ket{0},$$ (3)

with $\ket{0}$ the vacuum and $E$ the energy. Importantly, the poles of Eq. (3) are the energies of the zero quasi-momentum biexcitons, which we show in Fig. 3(g). Here, the moiré leads to gaps opening in the two-body (zero quasi-momentum) scattering continuum. We find that the higher energy biexcitons can exist in these gaps, consistent with our experimental observations.

To directly compare the theoretical model to the exciton shifts measured in experiment, we employ a polaron inspired model in which a single $\uparrow$ exciton is dressed by finite quasi-momentum excitations out of a coherent state of $\downarrow$ excitons generated by the pump. Employing a classical description of the pump light, we find the self-energy of the $\uparrow$ exciton within the so-called ladder approximation

$$\Sigma_{\lambda_{\uparrow}\lambda_{\uparrow}^{\prime}}(E;\omega_{L})=\sum_{\lambda_{\downarrow}\lambda_{\downarrow}^{\prime}}\beta_{\lambda_{\downarrow}}T^{\lambda_{\uparrow}\lambda_{\downarrow};\mathbf{0}}_{\lambda_{\uparrow}^{\prime}\lambda^{\prime}_{\downarrow};\mathbf{0}}(E+\omega_{L}+2i\gamma)\beta^{*}_{\lambda_{\downarrow}^{\prime}},$$ (4)

where the magnitude of $|\beta_{\lambda}|\propto\sqrt{n_{\lambda}}$ (with $n_{\lambda}$ the density of $\downarrow$ excitons in the $\lambda$ band; see Appendix H.4 for details) and $\gamma$ is the inverse lifetime of the excitons. Here, $\Omega$ is the light-matter coupling and we observe that the pump light frequency $\omega_{L}$ effectively leads to off-shell scattering of the excitons. The shift measured in the MX₁ ($\lambda=0$) and MX₂ ($\lambda=3$) excitons is given by $\Delta_{\lambda}(\delta_{1})\simeq\Sigma_{\lambda\lambda}(E_{\mathbf{0},\lambda};\omega_{L}=E_{\mathbf{0},0}-\delta_{1})$, where we have introduced the detuning $\delta_{1}$ as defined in experiment.

Figure 3(f) shows the ac Stark shifts predicted by our model. It is necessary to normalize the experimental light shifts by the pump laser power $l_{\rm pump}$. Consequently, we normalize the theoretical light shifts by $\hbar^{2}\Omega^{2}\propto l_{\rm pump}$ for comparison. Remarkably, our findings are in good qualitative agreement with the experimental observations. The theoretical calculation also suggests that for the measured experimental detunings, only two biexcitons lead to a noticeable sign change in the MX₁ and MX₂ energies.

III Attractive polaron light shift

Figure 4: Light shift of moiré polarons. (a) The peak reflection contrast as a function of gate voltage, showing that the oscillator strength of the AP linearly increases until $\nu=1$ and then decreases linearly until $\nu=2$. (b) Reflection spectrum $\Delta R_{\mathrm{AP}}(\tau)$ of the AP at $\nu=1$ as a function of the pump–probe delay $\tau$. (c) The observed AP blue shift is not driven by interactions as shown by its dependence on the pump detuning. (d) The measured light shift is remarkably constant as a function of $\nu$, up to a small dispersion $\sim 3\%$ (shown as a blue stripe) compatible with statistical fluctuations (shown as error bar, obtained for each points by repeating four time the measurement). The combination of linear $1/\delta_{\mathrm{AP}}$ dependence of the light shift and its independence of $\nu$ demonstrate the lack of interactions between APs or trions localized on different moiré sites.

Next, we focus on the attractive polaron (AP) resonance which emerges around $\approx 35\,$meV on the red side of MX₁. The corresponding trion binding energy is about $50\%$ larger than that observed in monolayer MoSe₂, suggesting that the electron Wannier orbital is strongly localized around the minimum of the moiré potential. The picture of localized moiré trion is also supported by the doping dependence of the AP oscillator strength, which follows the number of singly occupied moiré sites: It increases linearly until $\nu=1$ and then decreases linearly until $\nu=2$ (Fig. 4 a).

Figure 4 b shows the pump–probe measurement as a function of $\tau$ for $\delta_{\mathrm{AP}}=12$ meV, clearly showing a blue shift of the AP resonance with negligible alteration of the resonance for $\tau\geq 0.2$ ps. While this measurement is reminiscent of the result for MX₁, Fig. 4 c shows a striking difference in $\delta_{\mathrm{AP}}$ dependence; namely, the AP light shift is better fit using a $1/\delta_{\mathrm{AP}}$ dependence, typical of the ac Stark shift observed for an ensemble of non-interacting two-level emitters. We explain this observation by arguing that the AP resonance can be considered as stemming primarily from a collective excitation of trions localized at the M-M sites of the moiré lattice. In the absence of inter-site hopping, the excitation of a trion at a given site cannot depend on the existence of a trion on any other site, indicating that moiré trions are non-interacting excitations. The only contribution to the light shift in this limit will come from the ac-stark shift of each independent site, whose magnitude scales as $1/\delta_{\mathrm{AP}}$. Small but non-zero hybridization of the collective trion excitation with the bare exciton could give rise to a finite interaction strength and a deviation form the pure $1/\delta_{\mathrm{AP}}$ contribution to the light shift. We note that recent experiments on the same moiré structure yielded magnetization signatures consistent with unexpectedly weak inter-site hopping [7].

We find a confirmation of this explanation when we measure the electron density dependence of the light shift: Figure 4 d shows that varying $\nu$ from $0$ to $2$ results in negligible variation in the magnitude of the light shift, despite large variation in the pump-induced virtual AP population. This behavior contrasts with the strong electron density dependence of the AP light shift previously observed in a monolayer MoSe₂ [25], as well as that of MX₃, which we identify as an itinerant moiré exciton-polaron existing for filling $\nu>1$ (see Fig. 10 in the Appendix): in stark contrast to the AP, we observe a clear electron density dependence of MX₃. Figure 10, together with Figs. 4 c–d demonstrate the power of nonlinear pump–probe measurements in uncovering the strikingly different nature of localized moiré AP from that of itinerant exciton-polaron resonances.

IV Resonant excitation of the attractive polaron

Figure 5: Resonant excitation of moiré polarons. When the pump and probe lasers are cross-circularly polarized, they drive the AP transition at different moiré sites, hosting electrons with opposite spins (see the sketch in d). As a result, the pump has a negligible effect on the probe spectrum (a). A linearly polarized pump on the other hand, can generate an AP on all sites, leading to a blue shift and most prominently a reduction of the AP oscillator strength (b). For a sufficiently large intensity, all sites ends up in a statistical mixture of a single electron and a moiré trion and consequently the amplitude of the AP resonance is divided by a factor two around zero time delay (c). This phenomenon is observed at all fillings of the lattice (e). Here, the blue dots are the fitted AP amplitudes in absence of a pump pulse. The red squares show the amplitude after the pump pulse, for increasing pump power. We observe a saturation at half the reference amplitude (black dotted line).

Having identified the AP resonance as a collective excitation of independent two-level emitters where the ground and excited states correspond to an electron and a trion at a given site, respectively, we address the possibility of saturating the two-level emitters by resonantly driving the AP resonance. Since we are not interested in measuring the reflection of the strong resonant pump laser, we block it with a polarization filtering. We investigate both cross-circularly and cross-linearly polarized pump and probe configurations (Fig. 5 a, b). For cross-circular polarization, the pump generates APs in sites hosting an electron in a given valley, while the probe addresses the sites where electrons are in the opposite valley. Consequently, we observe no change in the probe reflection spectrum even for $\tau\simeq 0$.

In stark contrast, a linearly polarized pump pulse addresses all sites, leading to a blue shift of the AP resonance in conjunction with a drastic reduction of its oscillator strength. These two effects relax on a timescale of $\approx 5$ ps – much longer than the pump pulse duration $\tau_{\mathrm{pump}}\simeq 0.1$ ps. Such an unexpectedly long relaxation time may stem partly from a long radiative lifetime of trions and partly from generation of dark collective trion states [40, 41, 42, 43]. For short timescales $\tau\lesssim\tau_{\mathrm{pump}}$, we observe a strong saturation of the AP amplitude, so that the resonant probe reflection contrast reduces to $50\%$ of its value in the absence of the pump pulse (c). This observation indicates that the pump intensity is strong enough to ensure that each moiré site hosting a single electron is driven into a balanced mixture of its ground (single electron) and excited (trion) states. We observe these signatures for all electronic densities in the range $0\leq\nu\leq 2$ (Fig. 5 e).

The experiments detailed in Fig. 5 present a realization of a nonequilibrium Bose–Fermi mixture consisting of electrons in a flat moiré band and optically injected excitons [44, 45]. The choice of resonant excitation of the AP transition forces the mixture into a state that can be described as a high density moiré trion gas. Using a circularly polarized resonant pump laser and an external magnetic field to valley polarize electrons, it may be possible to create a trion at each moiré site, thereby realizing a solid-state analog of the Dicke model.

V Discussion

Our findings shed new light on the nature of moiré excitons and polarons, which remains a topic of active research [46, 47, 48, 49, 50, 5, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 32]. In particular, we show how non-linear spectroscopy unveil the itinerant or localized character of moiré optical excitations. This insight is crucial for the interpretation of experiments aimed at optical sensing of correlated electronic states [17, 18, 7] and the realization of degenerate Bose-Fermi mixtures in van der Waals heterostructures [44, 45]. Our experiments allow us to assess the extent of spatial overlap between different moiré excitons or attractive polarons. It is somewhat remarkable that this information is accessible to far-field optics given that the moiré length scale is about two orders of magnitude bellow the optical resolution.

The dressing of a quantum material with virtual optical excitations could be a promising route to engineer new phases of matter. The idea of using an intense laser pulse to modify material properties has already been demonstrated, but it often suffers from incoherent pumping and heating effects when an electronic polarization mode of the system is resonantly driven. Similar, albeit much less severe, problems in driven atomic systems can be alleviated using Rydberg dressing, where an off-resonant laser effects coherent hybridization of a ground (or long-lived low-energy) state and a Rydberg state: The atoms remain mostly in the ground state and spontaneous emission is strongly suppressed, but virtual excitations to the Rydberg state still ensure long-range interactions. The strategy we followed is somewhat analogous, albeit employed for a different purpose: The pump laser realizes an excitonic dressing of the moiré system, while keeping light absorption negligible thanks to a large enough detuning from electronic resonances. In turn, the probe pulse monitors the coherent response of the dressed moiré material, yieding information about the nature of the elementary excitations and their interactions.

VI Acknowledgement

We thank M. Hafezi, A. Srivastava, A. Christianen, A. G. Salvador and A. Müller for inspiring discussions. This work was supported by the Swiss National Science Foundation (SNSF) under Grant Number 200020$\_$207520. B. E. acknowledges funding from an ETH postdoc fellowship. H. S. A. acknowledges support from the Swiss Government Excellence Scholarship. T. U. acknowledges support from the Funai Overseas Scholarship.

The data are available at the ETH Research Collection [61].

VII Appendix

Appendix A Experimental setup and sample fabrication

The sample was measured in a dry cryostat (Attodry800, Attocube) at cryogenic temperaures of $\approx 5$K with free-space optical access and equipped with nanopositioners allowing displacement along the three axes. For the pump and probe, we used a mode-locked Ti:sapphire laser (Tsunami, Spectra-Physics), with a repetition rate of 76 MHz and pulse duration $\approx 100$ fs. The pulse is split along two paths for the pump and the probe. The bandwidth of the pump is reduced using a pulse shaper and its power is controlled using a motorized optical attenuator. We achieve a larger spectral width for the probe using a nonlinear fiber (femtowhite 800, NKT Photonics) which produces a quasi-continuum around the investigated resonances. The length of the probe optical path can be varied using a retroreflector on a motorized translation stage, enabling a fine tuning of the time delay between the two pulses. Both beams were focused on a diffraction-limited spot on the sample using an apochromatic microscope objective with NA $=0.8$ (LT-APO/VISIR/0.82, Attocube). The typical probe power is on the order of a few microwatts. The reflected light spectra were recorded using a Peltier-cooled CCD camera.

For the sample fabrication, few-layer graphite, $\approx 35$ nm hBN, monolayer MoSe₂ and WS₂ were mechanically exfoliated. The layers were assembled using the dry-transfer technique with a poly(bisphenol A carbonate) film on a polydimethylsiloxane (PDMS) stamp and deposited on a 285 nm Si/SiO₂ substrate [62]. The crystal alignement of the TMDs was determined prior to the stacking measuring the generation of second-harmonic light as a function of the polarization of an incoming infrared laser pulse. They were then stacked with a negligible twist. The graphene top and bottom gates and TMDs were contacted using gold electrodes deposited using optical lithography and electron beam deposition.

Appendix B Data analysis

In order to extract the light shift amplitude, we first fit the reflection spectrum for various time delays and extract the resonance position(s). The result is shown in Fig. 6 a  in the case of the MX₁ resonance. We then perform a Gaussian fit of the measured line shift as a function of the pump–probe delay in order to exctract the light shift amplitude at zero time delay. To determine the dependence of the light shift on the pump detuning, we perform for each pump wavelength a measurement at various pump powers as shown in Fig. 6 b. In this way, we can use relatively low power for a near resonant excitation and larger power at larger detunings, always making sure that we are in a regime of linear scaling with power. In the figures of the main text, we plot the slope $\Delta/I_{\rm pump}$ computed from a linear fit at each pump detuning $\delta$.

Figure 6: Data analysis. (a) Light shift as a function of the pump–probe delay for MX₁ for a pump power of $250\text{\,}\mathrm{\SIUnitSymbolMicro W}$ and detuning of 11 meV, showing a clear blue shift at zero time delay. The fit corresponds to a Gaussian envelope. (b) The dependence of the light shift of MX₁ on pump power at different pump detunings $\delta_{1}$ exhibits a linear dependence.

Appendix C Fit of the wavelength dependence

We discuss here our analysis of the wavelenght dependance of the light shift. Taking into account a single resonance, the light shift $\Delta$ can be expanded in a series of $1/\delta$ using perturbation theory [21]. The first order $1/\delta$ term comes from light-matter interaction, while exciton-exciton interactions contribute to higher order term, in $1/\delta^{2}$ or $1/(\delta_{i}\delta_{j})$ when several modes contributes. In principle, a fine analysis of the detuning dependence of the light shift would enable to deconvolve the various contributions. In practice, this analysis can be challenging due to the finite range of detuning in which we can take reliable data. Indeed, for all resonances, we had to restrict the data to a window of $\approx[20,80]$ meV. Going closer to resonance we face two issues: First, incoherent effects become more prominent as pump photons carry enough energy to generate a real population of excitons and hence are more likely to be absorbed. Second, the perturbative expansion of the light shift becomes inconsistent when $\Delta\sim\delta$. Conversely, going further away from resonance the signal reduces (given the available laser power) and the extraction of the light shift becomes unreliable. Furthermore, in that regime, the light–matter terms $\propto 1/\delta$ becomes stronger compared to the interaction term $\propto 1/\delta^{2}$ that we are interested in.

C.1 Charge neutrality

Figure 7: Pump wavelength dependence. From a fit we infer the origin of the light shift of MX₁ (a), MX₂ exciton (b) and of the attractive polaron at $\nu=1$ c. At charge neutrality (a and b), we are not able to deconvolve the contribution of the $u_{1,2}^{\prime}$ and $k_{2}^{\prime}$ fitting parameters (see Eq. [8]). We thus perform two fits, setting $k_{2}^{\prime}$ to its upper and lower bounds, shown as solid blue and red line, respectively. For the MX₁ exciton (a), the light shift is dominated by the $2u_{1,1}n_{1}$ term (shown as dashed line) and barely sensitive to the MX₂-dependent terms, so both fit are nearly identical and we obtain a good estimate of $u_{1,1}$. For MX₂ (b), both fits are in fair agreement with the data, but the contribution of the $2u_{1,2}n_{1}$ term (dotted line) is significantly different, and we are thus only able to provide a rather broad confidence interval for $u_{1,2}$. For the AP, a fit $A_{\mathrm{AP}}/\delta_{1}+B_{\mathrm{AP}}/\delta_{\mathrm{AP}}^{2}$ (solid black line, c) suggests that the usual AC-Stark shift scaling as $A_{\mathrm{AP}}/\delta_{1}$ (dashed-dotted line) is occurring, instead of an interaction-driven shift scaling as $B_{\mathrm{AP}}/\delta_{\mathrm{AP}}^{2}$ (dotted line).

With the relatively small detuning we use (compared to the exciton Rydberg energy), the exciton-exciton interaction is expected to be the dominant contribution to the light shift [25]. This is indeed confirmed by our data which are not compatible with a $1/\delta$ law. Let us thus focus on the interaction induced light shift, which we generalize to our moiré system with two bright exciton modes. As shown in App. H.6, the light shifts of MX₁ and MX₂ read

$$\Delta_{1} =2u_{1,1}n_{1}+2u_{1,2}n_{2}+4k_{1}\sqrt{n_{1}n_{2}}\,,$$ (5)

$$\Delta_{2} =2u_{2,2}n_{2}+2u_{1,2}n_{1}+4k_{2}\sqrt{n_{1}n_{2}}\,.$$ (6)

with the densities scaling as $n_{i}\propto|\Omega_{i}\phi_{i}(0)|^{2}/\delta_{i}^{2}$. The ratio $r=|\Omega_{1}\phi_{1}(0)|^{2}/|\Omega_{2}\phi_{2}(0)|^{2}=\gamma_{\rm rad,1}/\gamma_{\rm rad,2}\approx 0.48$ can be obtained from a fit of the reflection contrast (see App. D). Using this independent estimate and performing a joint fit of $\Delta_{1,2}$ we reduce the number of fitting parameters to five, $A$, $u_{2,2}^{\prime}$, $u_{1,2}^{\prime}$, $k_{1}^{\prime}$ and $k_{2}^{\prime}$, where $u_{i,j}^{\prime}=u_{ij}/u_{1,1}$, $k_{i}^{\prime}=k_{i}/u_{1,1}$ and

$$\Delta_{1} =A\left(\frac{1}{\delta_{1}^{2}}+\frac{r^{2}u_{1,2}^{\prime}}{\delta_{2}^{2}}+\frac{2rk_{1}^{\prime}}{\delta_{1}\delta_{2}}\right)\,,$$ (7)

$$\Delta_{2} =A\left(\frac{r^{2}u_{2,2}^{\prime}}{\delta_{2}^{2}}+\frac{u_{1,2}^{\prime}}{\delta_{1}^{2}}+\frac{2rk_{2}^{\prime 2,2}}{\delta_{1}\delta_{2}}\right)\,.$$ (8)

Furthermore, within a Born approximation $U_{ij}^{kl}\propto\int dr\phi_{i}\phi_{j}\phi_{k}^{*}\phi_{l}^{*}$ (see App H.6) and using the Cauchy-Schwarz inequality, we obtain the following bounds $0<u_{1,2}^{\prime}<[u_{2,2}^{\prime}]^{1/2}$ , $|k_{1}^{\prime}|<[u_{1,2}]^{1/2}$ , $|k_{2}^{\prime}|<[u_{2,2}^{\prime}u_{1,2}^{\prime}]^{1/2}$, which we enforce to improve the convergence of the fit.

For MX₁, see fig. 7 a, the full fit reveals the dominant contribution to be the MX₁-MX₁ interaction. We obtain the fitting parameter $A$, which we can compare to the value obtain from fitting the light shift of a monolayer exciton (on a monolayer MoSe₂ region of the same device). We estimate $u_{1,1}/u_{\rm ex}\approx 1.6\pm 0.2$, where $u_{\rm ex}$ is the monolayer exciton-exciton interaction (see App. E for an absolute calibration).

For MX₂, see fig. 7 b, the fit is unable to reliably disentangle the contribution of the second and third terms of 8, corresponding the sattering of an MX₁ exciton, without (1st term) or with (2nd term) a change of moiré band. Indeed, performing a fit with each of these two terms independently, we obtain in both case a reasonable agreement with the data over the full detuning range. We thus perform two fits, where the coupling $k_{2}^{\prime}=\pm[u_{2,2}^{\prime}u_{1,2}^{\prime}]^{1/2}$ staturates the Cauchy-Schwartz inequality. In this way we obtain the bounds $0.3\lessapprox u_{2,2}\lessapprox 0.8$.

C.2 Moiré attractive polaron

For the AP (c), the light shift is typically smaller and more noisy due to the weakness of the transition. As a result, it is more difficult to discriminate a potential $1/\delta$ and $1/\delta^{2}$ dependence. Nevertheless, the fit does suggest that the former is here the leading contribution, consistent with the picture of an ensemble of distinguishable and non-interacting two-level systems, as argued in the main text.

Appendix D Transfer matrix simulation

Resonance	$\nu$	$\hbar\gamma_{\rm rad}$ [meV]	$\hbar\gamma_{\text{non-rad}}$ [meV]
MX1	0	$1.0\pm 0.2$	$1.3\pm 0.3$
MX2	0	$0.48\pm 0.1$	$2.8\pm 0.4$
AP	1	$0.27\pm 0.1$	$2.7\pm 0.4$
MX${}_{2}^{\prime}$	1	$1.1\pm 0.3$	$2.9\pm 0.4$
MX3	2	$0.52\pm 0.1$	$2.2\pm 0.3$

Figure 8: Fit of the reflection spectrum using transfer matrix simulation. The blue dots are the data, the black solid line is a fit (grey area span by varying the fitting parameters within the confidence interval).

In order to infer the interaction strength of the various moiré excitons from their light shift, we need to estimate the exciton density that we generate, and hence the exciton oscillator strength. The latter can be obtained from a fit of the reflection spectrum. Such a fit needs to include the reflection of the electromagnetic field on the interfaces between the various dielectrics that make our van der Waals heterostructure. We do this using the transfer matrix method [63]. Two fitting parameters for the background are the hBN thickness $\approx 31\,$nm and $\approx 37\,$nm for the top and bottom layer and the hBN refractive index $n_{\rm hBN}\approx 2.15$ [64, 65]. Then, for each resonance, we have three additional fitting parameters, namely its energy, radiative and non-radiative decay rates. The results of this fit are shown in Fig. 8. We show here the bare reflection spectrum, obtained using a light source that is to a very good approximation spectrally flat in the energy range shown in the figure. We are unable to obtain a perfect fit of the background, using the hBN thicknesses and refractive index as free parameters. The discrepancy that we observe, especially on the edge of the spectrum could be due to chromatic aberrations (although we are using a microscope apochromatic objective to limit those). Nevertheless, in the center of the spectrum, we are able to reproduce our spectra very well at all fillings.

Appendix E Interaction strength

The radiative ($\gamma_{\rm r}$) and non-radiative ($\gamma_{\rm nr}$) decay rates can be used to extract the density of excitons $n_{\rm ex}$ induced by the pump. Using the optical Bloch equations within the adiabatic approximation we have [63]

$$n_{\rm ex}=\frac{2I_{\rm i}\gamma_{\rm r}}{(\gamma_{\rm r}+\gamma_{\rm nr})^{2}+\delta_{\rm ex}^{2}}\,,$$ (9)

where $I_{\rm i}$ is the photon flux. For reference, we first look at the exciton light shift $\Delta_{\rm ex}$ measured on a monolayer region of the sample, as a function of the exciton density, see Fig. 9. We observe a linear dependence $\Delta_{\rm ex}=U_{\rm ex}n_{\rm ex}$, from which we extract the exciton–exciton interaction strength $U_{\rm ex}\approx$0.06(3)\text{\,}\mathrm{\SIUnitSymbolMicro eV}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$$. The latter is compatible with previous measurements [66, 63, 40, 25], although significantly smaller than theoretical estimates, $\sim 3\,E_{\mathrm{ex}}a_{\mathrm{ex}}^{2}\sim$1\text{\,}\mathrm{\SIUnitSymbolMicro eV}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$$, where $E_{\mathrm{ex}}$ is the exciton binding energy and $a_{\mathrm{ex}}$ is the exciton Bohr radius [36, 37].

From the measurement of the light shift of the moiré excitons and polaron shown in Fig. 7 we obtain $U_{\text{MX1--MX1}}\approx 1.6\,U_{\rm ex}$ , $U_{\text{MX1--MX2}}\approx 0.7\,U_{\rm ex}$ and $U_{\text{MX3--MX3}}\approx 5.3\,U_{\rm ex}$. The enhancement of the interaction strength of MX₁ (charge neutrality) could stem from the partial confinement induced by the moiré potential. The larger enhancement for MX₃ reflects the polaronic nature of this resonance, as previously observed in a monolayer sample [25].

Figure 9: Light shift as a function of the exciton density. The blue dots are the data, obtained for various powers and detunings. The error bars come from the uncertainty of the exciton decay rates and on the power on the sample. The black solid line is a fit with the grey area span obtained by varying the fitting parameters within the confidence interval.

Appendix F Light shift of MX₃ in the regime of large doping

Figure 10: Light shift of MX₃. The resonance MX₃ shows distinctive kinks in its amplitude (a) and energy (b) at integer fillings of the moiré lattice with electrons, following the changes in the 2DES compressibility. The light shift of the moiré exciton MX₃ at various fillings (c), is well fitted by a line on top of a Gaussian function (black line). The latter captures the coherent response of the system, which corresponds to an interaction-induced blue shift for co-circularly polarized pump and probe, and possibly to a red shift stemming from the coupling to the biexciton states in a cross-circular configuration. The amplitudes of these two shifts are sensitive to the polaron dressing of the exciton, and consequently are extremal at integer filling of the moiré lattice (d).

Even though MX₃ is the dominant excitonic resonance for electron filling factors $\nu\geq 1.5$, its identification has remained unclear. To gain insight, we investigate the nonlinear response of the MX₃ resonance. Figure 10 a, b show the reflection amplitude and the energy of MX₃ for $1.5\leq\nu\leq 3.4$ in the absence of a pump laser: Consistent with earlier observations, we find that the resonance energy as well as the reflection strength, or equivalently the oscillator strength, of MX₃ exhibit local maxima at integer fillings $\nu=2$ and $\nu=3$. These features could be explained by partial suppression of dynamical dressing of excitons by electrons, when the two dimensional electron systems (2DES) is in an incompressible state [67, 17]. In contrast, when the electrons form a Fermi liquid ($\nu\neq 2,3$), the dynamical dressing of MX₃ is more effective and results in a red shift together with a reduction of the oscillator strength.

Figure 10 c shows the light shift as a function of pump–probe delay $\tau$ for four representative filling factors obtained for $\delta_{3}=80$ meV. We observe that the the light shift for $\tau\simeq 0$ indicates repulsive (attractive) interactions between same (opposite) valley MX₃ excitons generated by co- (cross-) circularly polarized pump–probe fields. While attractive interactions between opposite valley excitons has been reported before, it is surprising that the magnitude of the light shift is comparable in the two cases. The attractive interactions for the cross-polarized configuration may be explained through a near-resonant two-photon (pump+probe) excitation of the biexciton resonance at $\omega_{\mathrm{XX}}$. Verification of this hypothesis could be achieved by changing the pump detuning $\delta_{3}$ so as to probe both $\delta_{3}\leq\omega_{\mathrm{XX}}$ and $\delta_{3}\geq\omega_{\mathrm{XX}}$, since for the latter case, the biexciton-mediated interactions would become repulsive. As we had to choose $\delta_{3}<\omega_{\mathrm{XX}}$ to avoid strong background absorption, we could not verify the role of biexciton in the measured light shift.

In contrast to the light shift measurements in the charge-neutFral regime, we find that the pump pulse results in a MX₃ line shift that increases linearly with $\tau$ for $0.2\leq t\leq 3.0$ ps. Moreover, the linear shift at a given $\nu$ is identical for co- and cross-circularly polarized pump–probe configurations, but has a different sign for compressible ($\nu\neq 2,3$) and incompressible ($\nu=2,3$) electron states. We tentatively explain this feature by generation of free carriers by non-resonant absorption of pump-photons that change the electron density for timescales well exceeding the pump duration.

Since the MX₃ resonance energy has maxima (minima) for $\nu=2,3$ ($\nu=3/2,5/2$), any pump-induced change in electron density will result in a red (blue) shift of the resonance energy. While we do not understand why the red (blue) shift depends linearly on $\tau$ for $\tau>\tau_{\mathrm{pump}}$, we speculate that pump-induced charges are initially generated in high-energy bands and that they influence the nonlinear response only as they relax to the lowest-energy available moiré band.

We also observe in Fig. 10 c that the magnitude of the light shift for the co-circularly polarized pump–probe configuration is smaller for incompressible states. Plotting $\Delta_{3}$ as a function of $\nu$ (Fig. 10 d) shows that the blue shift is indeed minimal for $\nu=2,3$. This is at a first glance surprising given that the oscillator strength, and consequently the generated MX₃ population, is maximal for these incompressible states. However, it was recently shown that interactions between attractive exciton–polarons [25, 68, 40] mediated by their dressing cloud, are dramatically enhanced compared to those of bare excitons. Such an enhancement of interaction strength may overcome the reduction of the oscillator strength of MX₃ for $\nu\neq 2,3$. This tentative explanation suggests that the MX₃ mode may be identified as a second attractive polaron mode where the exciton is dressed by electrons in the upper moiré band. Last but not least, we find that the red shift of MX₃ in the cross-circularly polarized configuration is maximal when the electronic state is incompressible; we currently do not have an explanation for this observation.

Appendix G Light shift of MX${}^{\prime}_{2}$ at $\nu=1$

At a unity filling $\nu=1$ of the moiré potential, we observe two bright resonances, the AP which we discussed in detail in the main text, and MX${}_{1}^{\prime}$, emerging from MX₂, and which we now investigate. Depending on the pump and probe polarizations, we observed different behaviors. In co-circular polarization, we obtain the usual blue shift which we attribute to MX${}_{2}^{\prime}$–MX${}_{2}^{\prime}$ interactions. Contrary to other resonances, we cannot confirm this claim by an analysis of the detuning dependence. Indeed, we observed strong incoherent behavior for a pump blue detuned from the AP, and we therefore only explored the red detuned situation. Specifically, we explore the range $\delta_{\rm MX2^{\prime}}\approx[80,110]$ meV in which we observe no significant evolution of the light shift, as expected from a scaling as $1/\delta_{\rm MX_{2}^{\prime}}^{2}$, see Fig. 11 a. By contrast, in cross-circular polarization we observe a distinct redshift which diverges close to the AP resonance, in excellent agreement with a $1/\delta_{\rm AP}^{2}$ scaling and suggesting an attractive interaction between AP and MX${}_{2}^{\prime}$ in opposite valleys. We point out that a similar behavior was observed in a monolayer system, and tentatively attributed to the reduction of the phase space filling upon the generation of an AP, for an opposite valley exciton [25]. The pump power dependence (b) of the light shift shows an interesting behavior for a near-resonant excitation of the AP $\delta_{\rm AP}\approx 8$ meV. In that case, at high intensity, we expect a saturation of the AP density as the moiré potential is filled up, as described in the main text (although here, keeping a finite $\delta_{\mathrm{AP}}$ we are unable to fully saturate the transition). Indeed, we observe a sub-linear increase of the cross-polarized light shift attributed to AP–MX${}_{2}^{\prime}$ interactions. On the contrary, the co-polarized light shift which we attribute to MX${}_{2}^{\prime}$–MX${}_{2}^{\prime}$ interactions is linear in the pump power, as the density of MX${}_{2}^{\prime}$ remains far from saturation of the moiré lattice ($\delta_{{\rm MX}_{2}^{\prime}}\approx 80\,$meV).

Figure 11: Light shift of MX${}_{2}^{\prime}$ at $\nu=1$. (a) Detuning dependence of the light shift of the resonance MX${}^{\prime}_{2}$ at a filling $\nu=1$. Red dots correspond to cross-polarized pump and probe and show a red shift of the MX${}^{\prime}_{2}$ resonance which can be well fitted (red line) by a dependence $\propto 1/\delta_{\rm AP}^{2}$ stemming from MX${}^{\prime}_{2}$–AP interactions. Blue dots correspond to co-polarized pump and probe and show a blue shift of the MX${}^{\prime}_{2}$ resonance which can be well fitted (blue line) assuming MX${}_{2}^{\prime}$–MX${}^{\prime}_{2}$ interactions and a scaling $\propto 1/\delta_{\rm MX_{2}^{\prime}}^{2}$. (b) Lightshift of MX${}_{2}^{\prime}$ as a function of pump power in different polarizations. Cross-circular polarization shows a distinct red shift; the deviation from a linear dependence is due to saturation of the AP population at large powers. Co-linear polarization shows a smaller blue shift which is well in the linear regime. (c) MX${}_{2}^{\prime}$ shift vs. pump–probe time delay in co-circular (blue) for a pump power of $500\text{\,}\mathrm{\SIUnitSymbolMicro W}$ and cross-circular (red) polarizations for $50\text{\,}\mathrm{\SIUnitSymbolMicro W}$.

Appendix H Theoretical model

In this section, we primarily provide a detailed description of the theoretical model used to analyze the charge-neutral, cross-circular polarization data. The foremost goal of this model is to qualitatively capture the energy shifts observed in the MX₁ and MX₂ exciton resonances when subjected to the influence of a pump laser.

As discussed in the main text, the energy shifts can be interpreted as the dressing of the probe exciton (denoted by $\uparrow$) by virtual excitons ($\downarrow$) created by the pump. To elucidate this mechanism, we begin by examining the scattering processes involving two distinguishable excitons. Specifically, we calculate the exciton-exciton $T$ matrix in the presence of a moiré potential, which is inspired by a recent treatment of two interacting atoms in a 2D square optical lattice [39]. This approach allows us to capture the essential physics of the interaction under the influence of the pump. The dressing effect is then incorporated by applying non-self-consistent $T$ matrix theory.

We conclude the section with a brief discussion on a model for the co-circular polarization data. For ease of notation we set the reduced Planck’s constant and the system area and to unity ($\hbar=\mathcal{A}=1$).

H.1 Model for interactions of two distinguishable excitons

We begin by modelling the interactions between two distinguishable rigid (1s) excitons in the presence of moiré. The Hamiltonian we consider consists of four terms

$$H=H_{X\uparrow}+H_{X\downarrow}+H_{d}+V.$$ (10)

The first two terms represent the exciton Hamiltonians for each valley ($\sigma={\uparrow,\downarrow}$):

$$H_{X\sigma}=\sum_{\mathbf{K}}\epsilon_{\mathbf{K},X}\hat{X}^{\dagger}_{\mathbf{K}\sigma}\hat{X}_{\mathbf{K}\sigma}+\sum_{\mathbf{K}\mathbf{Q}}\tilde{V}_{X}(\mathbf{Q})\hat{X}^{\dagger}_{\mathbf{K}+\mathbf{Q}\sigma}\hat{X}_{\mathbf{K}\sigma}.$$ (11)

Here, $\hat{X}_{\mathbf{K}\sigma}$ annihilates an exciton of type $\sigma$ with energy $\epsilon_{\mathbf{K},X}=\epsilon_{X}+|\mathbf{K}|^{2}/2m_{X}$, where $m_{X}$ is the exciton mass and $\epsilon_{X}$ is the 1s exciton energy. The term $\tilde{V}_{X}$ is the Fourier transform of the exciton moiré potential, which is approximated in real space by [56]

$$V(\mathbf{x})=\sum_{j=1}^{6}V_{j}e^{i\mathbf{G}_{j}\cdot\mathbf{x}},$$ (12)

where $\mathbf{G}_{j}$ are the first six reciprocal lattice vectors. The three-fold symmetry and the realness of the periodic potential imply that $V_{1}=V_{3}=V_{5}$, $V_{2}=V_{4}=V_{6}$, and $V_{1}=V_{4}^{*}$, parametrized by $V_{1}=Ve^{i\psi}$, with $V$ determining the potential depth and $\psi$ its shape.

The third term in the Hamiltonian describes the closed-channel molecule, which mediates interactions between excitons:

$$H_{d}=\sum_{\mathbf{K}}(\epsilon_{\mathbf{K},d}+\delta_{\rm cc})\hat{d}^{\dagger}_{\mathbf{K}}\hat{d}_{\mathbf{K}}+\sum_{\mathbf{K}\mathbf{Q}}2\tilde{V}_{X}(\mathbf{Q})\hat{d}^{\dagger}_{\mathbf{K}+\mathbf{Q}}\hat{d}_{\mathbf{K}},$$ (13)

where $\hat{d}_{\mathbf{K}}$ annihilates a closed-channel molecule with energy $\epsilon_{\mathbf{K},d}=2\epsilon_{X}+|\mathbf{K}|^{2}/2M$ ($M=2m_{X}$) and detuning $\delta_{\rm cc}$. The closed-channel experiences both exciton potentials, hence the factor of two in the moiré potential.

The interactions are mediated according to

$$V=g\sum_{\mathbf{K}\mathbf{Q}}\chi(\mathbf{K})(\hat{d}^{\dagger}_{\mathbf{Q}}\hat{X}_{\mathbf{Q}-\mathbf{K},\uparrow}\hat{X}_{\mathbf{K},\downarrow}+\rm{h.c.}),$$ (14)

where $\chi(\mathbf{K})$ regularizes the ultraviolet (UV) divergence. Throughout this section we will use $\chi(\mathbf{K})=\Theta(\Lambda-|\mathbf{K}|)$, where $\Lambda$ is the UV cutoff. The bare coupling $g$ and detuning $\delta_{\rm cc}$ are renormalized as [69]

$$\frac{\delta_{\rm cc}}{g^{2}}=\sum_{\mathbf{K}}\chi(\mathbf{K})\frac{1}{E_{\rm BX}+2\epsilon_{\mathbf{K},X}},$$ (15)

with $E_{\rm BX}$ being the energy of the biexciton without moiré.

The two-channel model introduced above is often employed in the study of ultra-cold atoms. There it is understood as being inspired by the underlying microscopics, whereby the interactions between atoms are mediated by a molecular state (i.e., a Feshbach resonance). However, when taking the so-called single-channel model limit ($\delta_{\rm cc},g\to\infty$), the two-channel model is entirely equivalent to contact interactions with coupling constant $U=-g^{2}/\delta_{\rm cc}$, which connects the Hamiltonian presented here to the Hamiltonian in Eq. (II.3). Throughout this work, we will exclusively work in the single-channel model limit since, for our purposes, the two-channel model is only used as a tool to simplify the calculation of the exciton-exciton $T$ matrix.

H.2 Single-particle Hamiltonians

By Bloch’s theorem we can diagonalize the single-particle Hamiltonians in Eq. (10). For an exciton of type $\sigma$, we have

$$H_{X,\sigma}\ket{\mathbf{k},\lambda,\sigma}=E^{\sigma}_{\mathbf{k},\lambda}\ket{\mathbf{k},\lambda,\sigma}$$ (16)

and for the closed-channel molecule,

$$H_{d}\ket{\mathbf{k},\lambda,d}=(E^{d}_{\mathbf{k},\lambda}+\delta_{\rm cc})\ket{\mathbf{k},\lambda,d},$$ (17)

where $\lambda$ is the band index and $\mathbf{k}$ is the quasi-momentum. Here and throughout, quasi-momentum is denoted with lowercase letters, while real momentum is denoted with capital letters.

The Bloch states $\ket{\mathbf{k},\lambda,\sigma}$ and $\ket{\mathbf{k},\lambda,d}$ can be expanded according to

$$\ket{\mathbf{k},\lambda,\sigma} =\sum_{m,n}\phi^{(\mathbf{k},\lambda,\sigma)}_{(m,n)}\hat{X}^{\dagger}_{\mathbf{k}+\mathbf{G}_{mn},\sigma}\ket{0}$$ (18)

$$\ket{\mathbf{k},\lambda,d} =\sum_{m,n}\phi^{(\mathbf{k},\lambda,d)}_{(m,n)}\hat{d}^{\dagger}_{\mathbf{k}+\mathbf{G}_{mn}}\ket{0}$$ (19)

where $\mathbf{G}_{mn}=m\hat{G}_{0}+n\hat{G}_{1}$, with $\hat{G}_{0}$ and $\hat{G}_{1}$ the moiré reciprocal lattice basis vectors and $m$, $n$ integers. Here the expansion coefficients satisfy

$$E^{\sigma}_{\mathbf{k},\lambda}\phi^{(\mathbf{k},\lambda,\sigma)}_{\mathbf{m}} =\frac{|\mathbf{k}+\mathbf{G}_{\mathbf{m}}|^{2}}{2m_{X}}\phi^{(\mathbf{k},\lambda,\sigma)}_{\mathbf{m}}+\sum_{\mathbf{m}^{\prime}}\tilde{V}(\mathbf{G}_{\mathbf{m}^{\prime}})\phi^{(\mathbf{k},\lambda,\sigma)}_{\mathbf{m}+\mathbf{m}^{\prime}}.$$ (20)

where we have set $\mathbf{m}=(m,n)$. The equations for the dressed dimmer are almost identical, with $\sigma\to d$, $m_{X}\to M$ and $\tilde{V}\to 2\tilde{V}$.

Table 1 summarises the relevant parameters for our system as determined by large-scale Density Functional Theory (DFT) calculations [7]. Here, we also include the exciton moiré parameters which are determined by assuming a tightly bound electron-hole pair such that $V_{X}e^{i\psi_{X}}\equiv V_{e}e^{i\psi_{e}}+V_{h}e^{i\psi_{h}}$, where $V_{i}$ and $\psi_{i}$ are the moiré parameters for the exciton ($i=X$), electron ($i=e$) and hole ($i=h$).

Particle	$V$ [$\rm meV$]	$\psi$	$m$ [$m_{e}$]
Electron	$-6.3$	$0^{\circ}$	$0.45$
Hole	$1.9$	$59^{\circ}$	$0.55$
Exciton	$5.6$	$163^{\circ}$	$1$

Table 1: Parameters for electron, hole, and exciton from density functional theory calculations [7]. The mass of the particles is given in units of the bare electron mass $m_{e}$. The assumed moiré length is $a_{M}=8.2$ nm.

Using the DFT moiré parameters, we compare the theoretical exciton spectrum to experimental data. In particular, Fig. 12 shows the spectrum at charge neutrality, which qualitatively matches the experimental observations. The lowest two peaks in the spectrum correspond to MX₁ and MX₂. The exciton spectrum is given by

$$A(E)=-\frac{1}{\pi}\sum_{\lambda}\imaginary\left[\frac{|\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\uparrow)}|^{2}}{E-E_{\mathbf{0},\lambda}+i\gamma}\right],$$ (21)

which is proportional to the reflectance measured in experiment. Here, $\gamma=1$ meV is the inverse lifetime of the exciton, which we set to roughly match experiment. It is noteworthy that the DFT moiré parameters underestimate the separation between MX₁ and MX₂ by approximately 10 meV. We remind the reader, as shown in Fig. 3(g), there exists dark excitonic states between MX₁ and MX₂ (and above), which can play a significant role in the biexciton spectrum.

Figure 12: Exciton spectral function at zero doping using the DFT parameters in Table 1. In agreement with experiment, the theory predicts two dominant excitonic resonances, which we identify as MX₁ and MX₂. We have included a broadening from the finite exciton lifetime of $\gamma=1$ meV.

To conclude this section, we use the Bloch states introduced above to recast the two-channel interactions. Noting that the interactions will conserve quasi-momentum, the relevant matrix elements of the interactions are

$$\mathcal{V}^{\mathbf{q}\lambda_{d}}_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}}\equiv\bra{\mathbf{q},\lambda_{d},d}V\ket{\mathbf{q}-\mathbf{k},\lambda_{\uparrow},\uparrow;\mathbf{k},\lambda_{\downarrow},\downarrow}/g^{2}$$

$${}\quad=\sum_{\mathbf{N_{c}}\mathbf{N_{r}}}\chi(\mathbf{k}+\mathbf{N_{c}}/2-\mathbf{N_{r}})\phi^{(\mathbf{q},\lambda_{d},d)*}_{\mathbf{N_{c}}}\phi^{(\mathbf{q}-\mathbf{k},\lambda_{\uparrow},\uparrow)}_{\mathbf{N_{c}}/2+\mathbf{N_{r}}}\phi^{(\mathbf{k},\lambda_{\downarrow},\downarrow)}_{\mathbf{N_{c}}/2-\mathbf{N_{r}}}.$$ (22)

We point out that since the UV cutoff is large, we can take $\chi(\mathbf{k}+\mathbf{N_{c}}/2-\mathbf{N_{r}})\to\chi(\mathbf{N_{r}})$, where $\mathbf{N_{r}}$ can be understood as the relative reciprocal lattice vector of the two excitons.

H.3 T Matrix

To study the interactions we calculate the $T$ matrix, which provides an exact solution to the full two-body problem. To begin we introduce the free exciton and closed-channel Green’s functions:

$$\hat{G}^{(0)}(E) =\frac{1}{E-\hat{H}_{X\uparrow}-\hat{H}_{X,\downarrow}}$$ (23)

$$\hat{D}^{(0)}(E) =\frac{1}{E-\hat{H}_{d}}.$$ (24)

The two-body $T$ matrix is given by the infinite series

$$T =\hat{V}\hat{D}^{(0)}\hat{V}+\hat{V}\hat{D}^{(0)}\hat{V}\hat{G}^{(0)}\hat{V}\hat{D}^{(0)}\hat{V}+\dots$$

$$=\hat{V}\left(\hat{D}^{(0)}+\hat{D}^{(0)}\hat{V}\hat{G}^{(0)}\hat{V}\hat{D}^{(0)}+\dots\right)\hat{V}$$

$$=\hat{V}\hat{D}\hat{V}$$ (25)

where we have suppressed the energy dependencies for brevity and $\hat{D}(E)$ is the closed-channel Green’s function. Thus, by finding the closed-channel Green’s function, we can immediately calculate the $T$ matrix. This approach simplifies the calculation and provides easier access to the biexciton energies, which are both the poles of the $T$ matrix and $\hat{D}(E)$.

The closed-channel Green’s function is given by

$$\hat{D}(E) =\hat{D}^{(0)}+\hat{D}^{(0)}\hat{V}\hat{G}^{(0)}\hat{V}\hat{D}^{(0)}+\dots$$ (26a)

$$=\frac{1}{[\hat{D}^{(0)}(E)]^{-1}-g^{2}\hat{\Pi}(E)},$$ (26b)

where we have introduce the polarization bubble, which has matrix elements

$$\Pi^{\mathbf{q}}_{\lambda\lambda^{\prime}}= \sum_{\mathbf{k},\lambda_{\uparrow},\lambda_{\downarrow}}\mathcal{V}^{\mathbf{q},\lambda}_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}}\frac{1}{E-E_{\mathbf{k},\mathbf{q};\lambda_{\uparrow},\lambda_{\downarrow}}}\mathcal{V}^{\mathbf{q},\lambda^{\prime}*}_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}},$$ (27)

where $E_{\mathbf{k},\mathbf{q};\lambda_{\uparrow},\lambda_{\downarrow}}\equiv E^{\uparrow}_{\mathbf{q}-\mathbf{k},\lambda_{\uparrow}}+E^{\downarrow}_{\mathbf{k},\lambda_{\downarrow}}$.

We note that the polarization bubble calculation is numerically demanding due to the need to contract large tensors. In order to optimize the calculation, we reduce the size of interaction tensor $\mathcal{V}^{\mathbf{q},\lambda}_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}}$ by eliminating all elements that are smaller than a specified tolerance. We have checked that this tolerance does not affect the final calculation.

H.4 Effect of the pump laser

To incorporate the effects of the pump laser, which generates $\downarrow$ excitons in our notation, we consider the following modification to the Hamiltonian $\hat{H}_{X\downarrow}$ [21]

$$H_{X\downarrow} =\sum_{\mathbf{k},\lambda}E^{\downarrow}_{\mathbf{k},\lambda}\hat{X}^{\dagger}_{\mathbf{k},\lambda,\downarrow}\hat{X}_{\mathbf{k},\lambda,\downarrow}+\Omega{}^{*}e^{i\omega_{L}t}\sum_{\lambda}\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)*}\hat{X}^{\dagger}_{\mathbf{0},\lambda,\downarrow}$$

$${}\qquad+\Omega e^{-i\omega_{L}t}\sum_{\lambda}\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)}\hat{X}_{\mathbf{0},\lambda,\downarrow}.$$ (28)

Here, $\Omega$ represents the light-matter interaction strength; we assume that the light couples only to the zero-momentum exciton and use the rotating wave approximation ($|\epsilon_{X}-\omega_{L}|\ll\epsilon_{X}+\omega_{L}$). To remove the time-dependence, we apply the unitary transformation

$$\hat{U}(t)=\exp(-i\omega_{L}\sum_{\mathbf{k},\lambda}\hat{X}^{\dagger}_{\mathbf{k},\lambda,\downarrow}\hat{X}_{\mathbf{k},\lambda,\downarrow})$$ (29)

to move into the rotating frame. This gives us the Hamiltonian

$$\tilde{H}_{X\downarrow} =\sum_{\mathbf{k},\lambda}\left(E^{\downarrow}_{\mathbf{k},\lambda}-\omega_{L}\right)\hat{X}^{\dagger}_{\mathbf{k},\lambda,\downarrow}\hat{X}_{\mathbf{k},\lambda,\downarrow}$$

$$+\Omega{}^{*}\sum_{\lambda}\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)*}\hat{X}^{\dagger}_{\mathbf{0},\lambda,\downarrow}+\Omega\sum_{\lambda}\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)}\hat{X}_{\mathbf{0},\lambda,\downarrow}.$$

This Hamiltonian can be diagonalized using the multi-mode shift operator,

$$\hat{\underline{D}}(\beta)=\exp[\sum_{\lambda}\left(\beta_{\lambda}\hat{X}^{\dagger}_{\mathbf{0},\lambda,\downarrow}-\beta_{\lambda}^{*}\hat{X}_{\mathbf{0},\lambda,\downarrow}\right)],$$ (30)

which has the property that $\hat{\underline{D}}(\beta)\hat{X}_{\mathbf{0},\lambda,\downarrow}\hat{\underline{D}}^{\dagger}(\beta)=\hat{X}_{\mathbf{0},\lambda,\downarrow}-\beta_{\lambda}$. The diagonalized Hamiltonian is

$$\tilde{H}_{X\downarrow} =\sum_{\lambda}\left(E^{\downarrow}_{\mathbf{0},\lambda}-\omega_{L}\right)\left(\hat{X}^{\dagger}_{\mathbf{0},\lambda,\downarrow}-\beta^{*}_{\lambda}\right)\left(\hat{X}_{\mathbf{0},\lambda,\downarrow}-\beta_{\lambda}\right)$$

$${}\quad+\sum_{\mathbf{k}\neq\mathbf{0},\lambda}\left(E^{\downarrow}_{\mathbf{k},\lambda}-\omega_{L}\right)\hat{X}^{\dagger}_{\mathbf{k},\lambda,\downarrow}\hat{X}_{\mathbf{k},\lambda,\downarrow}+E_{0},$$ (31)

where $\beta_{\lambda}=-\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)*}\Omega^{*}/(E^{\downarrow}_{\mathbf{0},\lambda}-\omega_{L})$ and

$$E_{0}=\sum_{\lambda}\left[(E^{\downarrow}_{\mathbf{0},\lambda}-\omega_{L})|\beta_{\lambda}|^{2}+2\real\left(\Omega\,\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)}\right)\right]$$ (32)

is the ground state energy. Enforcing that all $\left(\hat{X}_{\mathbf{k},\lambda,\downarrow}-\delta_{\mathbf{k},\mathbf{0}}\beta_{\lambda}\right)$ annihilate the ground state ($\ket{\Phi}$), we find $\ket{\Phi}=\hat{\underline{D}}(\beta)\ket{0}$.

Through this analysis, we see that a classical treatment of the laser can be accounted for by substituting $\hat{X}_{\mathbf{0},\lambda,\downarrow}\to\hat{X}_{\mathbf{0},\lambda,\downarrow}-\beta_{\lambda}$ and $E^{\downarrow}_{\mathbf{k},\lambda}\to E^{\downarrow}_{\mathbf{k},\lambda}-\omega_{L}$. In principle, moving into the rotating frame causes the interaction in Eq. (14) to develop a time-dependence. However, we are free to cancel this time-dependence with a counter-rotating frame for the $\uparrow$ exciton terms, (i.e., $\downarrow\to\uparrow$ and $\omega_{L}\to-\omega_{L}$ in Eq. (29)). This will of course also shift the $\uparrow$ exciton spectrum by the constant $\omega_{L}$, but since we are only interested in physics measured relative to $\uparrow$ spectrum, this is not necessary to include.

H.5 Dressing probe exciton by pump excitons

We now consider how excitons generated by the pump (denoted by $\downarrow$) dress a single probe exciton ($\uparrow$). This is analogous to the theory of the Bose polaron, where we now effectively have multiple non-interacting Bose-Einstein condensates of $\downarrow$ excitons in various zero quasi-momentum bands. Within the framework of non-self-consistent T-matrix theory, the self energy of the probe exciton is given by

$$\Sigma_{\lambda_{\uparrow}\lambda_{\uparrow}^{\prime}}(E;\omega_{L})=\sum_{\lambda_{\downarrow}\lambda_{\downarrow}^{\prime}}\beta_{\lambda_{\downarrow}}T^{\lambda_{\uparrow}\lambda_{\downarrow};\mathbf{0}}_{\lambda_{\uparrow}^{\prime}\lambda^{\prime}_{\downarrow};\mathbf{0}}(E+\omega_{L}+2i\gamma)\beta^{*}_{\lambda_{\downarrow}^{\prime}},$$ (33)

where $T^{\lambda_{\uparrow}\lambda_{\downarrow};\mathbf{0}}_{\lambda_{\uparrow}^{\prime}\lambda^{\prime}_{\downarrow};\mathbf{0}}(E)\equiv\bra{\mathbf{0},\lambda_{\uparrow},\uparrow;\mathbf{0},\lambda_{\downarrow},\downarrow}\hat{T}(E)\ket{\mathbf{0},\lambda_{\uparrow}^{\prime},\uparrow;\mathbf{0},\lambda_{\downarrow}^{\prime},\downarrow}$ and both excitons contribute a term $i\gamma$ due to their finite lifetimes. Here, we explicitly include the laser frequency $\omega_{L}$ as a parameter in the self energy to emphasize its significance: we observe that the laser light shifts the scattering off-shell, similar to the effect discussed in the context of polariton-electron scattering [70].

Non-self-consistent $T$ matrix theory is a non-perturbative approach in that it sums over all ladder diagrams, capturing an infinite series of scattering events to account for strong interactions. Alternatively, the non-self-consistent $T$ matrix theory can be understood as entirely equivalent to a variational approach based on a generalization of the celebrated Chevy ansatz [71]:

$$\ket{\Psi} =\sum_{\lambda_{\uparrow}}\psi^{(0)}_{\lambda_{\uparrow}}\hat{X}^{\dagger}_{\mathbf{0},\lambda_{\uparrow},\uparrow}\ket{\Phi}+\sum_{\lambda}\psi_{\lambda}^{(0;d)}\hat{d}^{\dagger}_{\mathbf{0},\lambda}\ket{\Phi}$$

$${}\qquad\sum_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}}\psi^{(1)}_{\mathbf{k}\lambda_{\uparrow}\lambda_{\downarrow}}\hat{X}^{\dagger}_{\mathbf{k},\lambda_{\uparrow},\uparrow}\hat{X}^{\dagger}_{-\mathbf{k},\lambda_{\downarrow},\downarrow}\ket{\Phi}.$$ (34)

By deriving the variational equations, one can show that they can be rearranged into the form

$$E\psi^{(0)}_{\lambda_{\uparrow}}=E_{\mathbf{0},\lambda_{\uparrow}}\psi^{(0)}_{\lambda_{\uparrow}}+\sum_{\lambda_{\uparrow}^{\prime}}\Sigma_{\lambda_{\uparrow},\lambda_{\uparrow}^{\prime}}(E;\omega_{L})\psi^{(0)}_{\lambda_{\uparrow}^{\prime}},$$ (35)

which clearly demonstrates the equivalence of the two methods.

Since we only aim for qualitative accuracy, we approximate the shift in the exciton resonance by the self-energy,

$$\Delta_{\lambda}(\omega_{L})\simeq\Sigma_{\lambda\lambda}(E_{\mathbf{0},\lambda};\omega_{L}).$$ (36)

It is important to note that Eq. (36) captures the ac Stark effect, whereby the shift can change sign when the laser frequency $\omega_{L}$ is tuned into resonance with the biexciton bound state, i.e., $E_{\mathbf{0},\lambda}=\omega_{L}-E_{\rm{BX};\alpha}$, where $E_{\rm{BX};\alpha}$ denotes the energy of a zero quasi-momentum biexciton state. However, unlike the simpler case of a monolayer TMD, this sign change does not occur at all biexciton energies, as not all biexciton states will couple to the $\lambda$-th exciton. We can therefore expect in the experiment, which only observes shifts in the optically active MX₁ and MX₂, to not detect all biexciton states. This is discussed further in the main text.

H.6 Same valley interactions

To conclude our work we briefly consider interactions from excitons in the same valley, which do not support a bound state. Owing to this, we simply take the Born approximation of the $T$ matrix assuming repulsive contact interactions. Focusing on zero quasi-momentum, the interactions between the excitons are given by

$$\hat{\mathscr{V}}=\frac{g}{2}\sum U^{\lambda_{1}\lambda_{2}}_{\lambda_{3}\lambda_{4}}\hat{X}^{\dagger}_{\mathbf{0},\lambda_{1},\downarrow}\hat{X}^{\dagger}_{\mathbf{0},\lambda_{2},\downarrow}\hat{X}_{\mathbf{0},\lambda_{3},\downarrow}\hat{X}_{\mathbf{0},\lambda_{4},\downarrow}.$$ (37)

Here $g$ is the strength of the repulsion, and we have introduce the overlap integral

$$U^{\lambda_{1}\lambda_{2}}_{\lambda_{3}\lambda_{4}}=\int d^{2}r\,\varphi_{\lambda_{1},\downarrow}^{*}(\mathbf{r})\varphi_{\lambda_{2},\downarrow}^{*}(\mathbf{r})\varphi_{\lambda_{3},\downarrow}(\mathbf{r})\varphi_{\lambda_{4},\downarrow}(\mathbf{r})$$ (38)

where $\varphi_{\lambda,\downarrow}$ is the Bloch wavefunction in real space at zero quasi-momentum with band index $\lambda$ and valley index $\downarrow$. In order to derive the mean-field Hamiltonian given in the main text, we restrict ourselves to the optically bright excitonic states MX₁ ($\lambda=0$) and MX₂ ($\lambda=3$); for further details on these states see the discussion below Eq. (II.3). We furthermore use the fact that $\beta_{\lambda}=\kappa_{\lambda}\sqrt{n_{\lambda}}$ (with $\kappa_{\lambda}=e^{i\rm{arg}(-\phi_{\mathbf{0}}^{(\mathbf{0},\lambda,\downarrow)}{}^{*}\Omega^{*})}$) and $\bra{\Phi}\hat{X}^{\dagger}_{\mathbf{0},\lambda,\downarrow}\hat{X}_{\mathbf{0},\lambda^{\prime},\downarrow}\ket{\Phi}=\beta_{\lambda}^{*}\beta_{\lambda^{\prime}}$. Decomposing the interaction in Eq. (37) using mean-field Hartree-Fock theory then yields Eq. (1). We point out that for contact interactions with bosons, the Hartree and Fock terms are identical in magnitude and sign. The relevant overlap integrals are

$$u_{11}/g =U^{0,0}_{0,0}=2.2$$ (39a)

$$u_{22}/g =U^{3,3}_{3,3}=1.9$$ (39b)

$$u_{12}/g =U^{0,3}_{0,3}=1.3$$ (39c)

$$k_{1}/g =\kappa_{0}^{3}\kappa_{3}U^{0,0}_{0,3}=-1.2$$ (39d)

$$k_{2}/g =\kappa_{0}\kappa_{3}^{3}U^{3,3}_{3,0}=-0.4,$$ (39e)

where we have assumed that $\kappa$ are real. All other interaction terms can be derived from these five, since the overlap integral only depends only on the number of each index (e.g., $U^{1,1}_{1,0}=U^{1,1}_{0,1}=U^{1,0}_{1,1}=U^{0,1}_{1,1}$).