Interspecies exciton interactions lead to enhanced nonlinearity of dipolar excitons and polaritons in MoS2 homobilayers

Nonlinear interactions between excitons strongly coupled to light are key for accessing quantum many-body phenomena in polariton systems. Atomically-thin two-dimensional semiconductors provide an attractive platform for strong light-matter coupling owing to many controllable excitonic degrees of freedom. Among these, the recently emerged exciton hybridization opens access to unexplored excitonic species, with a promise of enhanced interactions. Here, we employ hybridized interlayer excitons (hIX) in bilayer MoS2 to achieve highly nonlinear excitonic and polaritonic effects. Such interlayer excitons possess an out-of-plane electric dipole as well as an unusually large oscillator strength allowing observation of dipolar polaritons (dipolaritons) in bilayers in optical microcavities. Compared to excitons and polaritons in MoS2 monolayers, both hIX and dipolaritons exhibit ≈ 8 times higher nonlinearity, which is further strongly enhanced when hIX and intralayer excitons, sharing the same valence band, are excited simultaneously. This provides access to an unusual nonlinear regime which we describe theoretically as a mixed effect of Pauli exclusion and exciton-exciton interactions enabled through charge tunnelling. The presented insight into many-body interactions provides new tools for accessing few-polariton quantum correlations.

when hIX and intralayer excitons, sharing the same valence band, are excited simultaneously.This gives rise to a highly nonlinear regime which we describe theoretically by introducing a concept of hole crowding.
An attractive approach for realization of dipolar excitons and polaritons is to employ the recently discovered exciton hybridization in MoS 2 bilayers [11,37].This  (4 K) showing three distinct absorption features at 1.937 eV, 2.004 eV and 2.113 eV for X A , hIX and hX B , respectively.The measured linewidths for X A , hIX, and hX B are 20, 23 and 64 meV, respectively.RC is calculated using the formula in the top-right corner of the graph.d, Sketch of the conduction and valence bands in two adjacent layers of MoS 2 , displaying the allowed optical transitions of A and B direct intralayer excitons (X A and X B ) and interlayer excitons (IX) for spin-up states (black lines) at the K point in the bilayer momentum space.IX hybridizes with X B through the hole tunnelling between the two layers (red dashed arrow).At the K' point of the bilayer Brillouin zone, the same configuration applies for the states with the opposite spins.e, RC spectra of excitons in BL MoS 2 detected in two circular polarizations in an out-of-plane magnetic field of 8 T at T =4 K. Zeeman shifts of opposite signs are observed for X A and hIX.The absorption peak of the charged intralayer exciton (X * A ) shows near unity circular polarization.approach allows realization of uniform samples suitable for the observation of macroscopic many-body phenomena [38].Interlayer excitons unique to bilayer MoS 2 possess a large oscillator strength, comparable to that of the intralayer exciton, arising from interlayer hybridization of valence band states, aided by a favourable orbital overlap and a relatively small spin-orbit splitting among semiconducting TMDs [11].
Such hybridized interlayer excitons (hIX) are highly tunable using out-of-plane electric field [12,13] and their valley degree of freedom persists up to room temperature [14].
Here we use hIXs in bilayer MoS 2 to realize highly nonlinear excitonic and dipolaritonic effects.We unravel a previously unexplored interaction regime involving intra-and interlayer excitons stemming from the fermionic nature of the charge carriers in a valence band shared between different excitonic species.This regime, accessible using broadband excitation resonant with both hIX and intralayer exciton transitions, provides strong (up to 10 times) enhancement of the exciton nonlinearity, already enhanced by up to 8 times in MoS 2 bilayers compared with monolayers.
We support our experimental findings with microscopic theory, analysing the excitonic many-body physics and the cross-interactions and introducing the nonlinear mechanisms of the hole crowding.
Our heterostructure samples consists of a MoS 2 bilayer (BL) sandwiched between hBN and placed on a distributed Bragg reflector (DBR).Fig. 1a shows a bright field microscope image of the encapsulated BL MoS 2 .A sketch of the side view of the device is displayed in Fig. 1b.The reflectance contrast (RC) spectrum of the studied MoS 2 bilayer, displayed in Fig. 1c, shows three peaks: the intralayer neutral excitons X A at at 1.937 eV (see Fig. 1d), hybridized interlayer exciton hIX at 2.004 eV and hybridized B-exciton at 2.113 eV.Due to the quantum tunnelling of holes, B-excitons hybridize with an interlayer exciton (IX) (Fig. 1d), which is a direct transition in the bilayer momentum space [11].The ratio of the integrated intensities of X A and hIX is 4.5.Based on these data, we estimate the electronhole separation d = 0.55 nm (see details in Supplementary Note S1) in agreement with previous studies [14].We further confirm the nature of the hIX states by placing the BL MoS 2 in magnetic field where the valley degeneracy is lifted (Fig. 1e).In agreement with recent studies [13,39], we measure a Zeeman splitting with an opposite sign and larger magnitude in hIX compared with X A (-3.5 versus 1.5 meV).
We study the strong coupling regime in a tunable planar microcavity (Fig. 2a) formed by a silver mirror and a planar DBR [7].RC scans as a function of the cavity mode detuning ∆ = E cav − E exc , where E cav and E exc are the cavity mode and the corresponding exciton energy, respectively, are shown in Fig. 2c,d.Characteristic anticrossings of the cavity mode with X A and hIX are observed, resulting in lower, middle and upper polariton branches (LPB, MPB, and UPB, respectively).The extracted Rabi splittings are Ω X A = 38 meV for X A and Ω hIX = 19 meV for hIX (Supplementary Note S2).Fig. 2d shows the RC spectra in the vicinity of the anticrossing with hIX, providing a more detailed view of the formation of the MPB and UPB.The intensity of the polariton peaks is relatively low for the states with a high exciton fraction at positive (negative) cavity detunings for the MPB (UPB).
As the Rabi splitting scales as a square root of the oscillator strength, the ratio Ω X A /Ω hIX = 2 is in a good agreement with the RC data for integrated intensities of X A and hIX.From the Rabi splitting ratio we can estimate the tunneling constant J leading to the exciton hybridization.The corresponding coefficient is J = 48 meV (see Supplementary Note S1 for details), matching the density functional theory predictions [11].In polarization-resolved cavity scans in an out-of-plane magnetic field (Fig. 2e,f), similarly to hIX behaviour, we observe opposite and larger Zeeman splitting for dipolaritons relative to the intralayer polaritons (see Supplementary Figure S4).Chiral dipolariton states are observed distinguished by their opposite circular polarization (Fig. 2f).
We investigate the nonlinear response of X A and hIX in the bare BL flake as a function of the laser power using both narrow band (NB, full-width at half maxi-mum, FWHM=28nm) and broad band (BB, FWHM=50 nm) pulsed excitation (see Methods).Our resonant pump-probe experiments have confirmed that the lifetimes of the hIX and X A states are considerably longer than the pulse duration of ≈ 150 fs (Supplementary Note S3).Measured RC spectra are shown in Fig. 3a,b for the NB and in Fig. 3c for BB excitation.In the NB case, the excitation was tuned to excite either X A or hIX independently, while in the BB case, both resonances were excited simultaneously.
As seen in Figs.3a,b both X A and hIX spectra behave similarly upon increasing the power of the NB excitation: a blueshift of several meV is observed, accompanied by the peak broadening and bleaching.For the BB excitation, however, a different nonlinear behaviour is observed as shown in Fig. 3c: the broadening and complete suppression of the hIX peak is observed at much lower powers, accompanied by a redshift.This is in contrast to X A , whose behaviour is similar under the two excitation regimes.
The resulting energy shifts, peak linewidths and intensities are shown in Fig. 3d,e as a function of the exciton density (see details in Supplementary Note S4 and S6).Fig. 3d quantifies the trends observed in Figs.3a,b showing for the BB excitation an abrupt bleaching of the hIX peak above the hIX density 5 × 10 3 µm −2 accompanied by a redshift of ≈ 4 meV and a 12 meV broadening.For the NB case, a similar decrease in peak intensity is observed only around 4×10 4 µm −2 , accompanied with a peak blueshift of ≈ 7 meV and a broadening exceeding 15 meV.In Fig. 3e, however, it is apparent that the observed behaviour under the two excitation regimes is similar for X A .A similar blueshift, broadening and saturation are observed at slightly higher densities compared to the hIX under the NB excitation (Supplementary Note S5).We also find that due to the increased excitonic Bohr radius, the onset of the nonlinear behaviour for X A in bilayers occurs at a lower exciton density than for X A in monolayers (Supplementary Note S7).
We develop a microscopic model to describe the contrasting phenomena under the NB and BB excitation.Under the NB excitation, either X A or hIX excitons are created as sketched in Fig. 4a.In this case, nonlinearity arises from Coulomb exciton-exciton interactions causing the blueshift and dephasing [40].For simplicity, in the main text we will use a Coulomb potential V Coul combining the exchange and direct terms further detailed in Supplementary Note S8.We confirm (see Supplementary Note S8) that for the intralayer exciton-exciton interaction (X A -X A ) the dominant nonlinear contribution comes from the Coulomb exchange processes, as in the monolayer case [40,41], while for the hIX-hIX scattering the dominant contribution is from the direct Coulomb (dipole-dipole) interaction terms [19].For both X A and hIX, the Coulomb interaction is repulsive, and thus leads to the experimentally observed blueshifts.We find that for the modest electron-hole separation d = 0.55 nm in the bilayer, V Coul is overall 2.3 times stronger for hIX compared with X A .
Analysing the shapes of the reflectance spectra in the NB case, we note that they depend on the rates of radiative (Γ R ) and non-radiative (Γ NR ) processes.The area under RC curves is described by the ratio Γ R /(Γ R + Γ NR ).This ratio changes under the increased excitation if the rates depend on the exciton densities.Specifically, we account for the scattering-induced non-radiative processes that microscopically scale as Γ NR ∝ |V Coul | 2 n, i.e. depend on the absolute value of the combined matrix elements for the Coulomb interactions and the exciton density n [40].This process allows reproducing the RC behaviour and bleaching at increasing pump intensity.
Moreover, it explains stronger nonlinearity for X A in bilayers compared to monolayers.Namely, the scattering scales with the exciton Bohr radius, V Coul ∝ α, which is larger in the bilayers due to the enhanced screening (Supplementary Note S8).
In the BB case, both X A and hIX excitons are generated simultaneously, and together with intraspecies scattering (X A -X A and hIX-hIX), interspecies scattering (X A -hIX) occurs, similarly to the direct-indirect exciton Coulomb scattering in double quantum wells [42].Since X A and hIX are formed by the holes from the same valence band (Fig. 4a), an additional contribution arises from the phase space filling, i.e. the commutation relations for the excitons (composite bosons) start to deviate from the ideal weak-density limit once more particles are created [43].For particles of the same flavour, the phase space filling enables nonlinear saturation effects in the strong coupling regime, similar to polariton saturation observed in [29].
However, in the presence of several exciton species, we reveal a distinct phase space filling mechanism which we term the hole crowding.Crucially, we observe that the commutator of the X A annihilation operator ( X) and hIX creation operator ( Î † ) is non-zero, [ X(p), Î † (q)] = − Bp,q .Here p, q are exciton momenta and Bp,q is an operator denoting the deviation from the ideal commuting case ( Bp,q = 0) of distinct bosons where holes do not compete for the valence band space.
This statistical property of modes that share a hole has profound consequences for the nonlinear response.Namely, the total energy is evaluated as an expectation value over a many-body state with both X A and hIX excitons, |N X , N hIX := where N X and N hIX particles are created from the ground state |Ω max .If the excitonic modes are independent, the contributions from X A and hIX simply add up.However, the hole coexistence in the valence band induces the excitonic interspecies scattering.The phase space filling combined with the Coulomb energy correction leads to a negative nonlinear energy contribution.This nonlinear term scales as ∆E hIX = −η √ n X n hIX , where η > 0 is a coefficient defined by the Coulomb energy and Bohr radii and n X,hIX are the exciton densities (see Supplementary Note S9).This nonlinearity also modifies the non-radiative processes leading to substantial broadening for the hIX states.
According to this analysis, the effect of the BB excitation should be most pronounced for hIX.In addition to the possible hIX-hIX scattering (similar to that occurring under the NB excitation), much stronger X A absorption leads to the phase space filling in the valence band.Such hole crowding introduces additional scattering channels for hIX and leads to its RC spectra bleaching at lower hIX exciton densities.On the other hand, as only relatively small hIX densities can be generated, both the NB and BB excitation cases should produce similar results for   We investigate nonlinear properties of dipolar polaritons in a monolithic (fixedlength) cavity created by a silver mirror on top of a PMMA spacer (245 nm thick) covering the hBN-encapsulated MoS 2 homobilayer placed on the DBR.The cavity mode energy can be tuned by varying the angle of observation (0 degrees corresponds to normal incidence).We use a microscopy setup optimized for Fourier-plane imag-ing, thus allowing simultaneous detection of reflectivity spectra in a range of angles as shown in Fig. 5(a) displaying the measured polariton dispersion.In this experiment, the cavity mode is tuned around hIX and only two polariton branches LPB and UPB are observed at low fluence of 0.6 µJ cm −2 with a characteristic Rabi splitting of 17.5 meV.In Fig. 5(b), at an increased fluence of 58.5 µJ cm −2 , only a weakly coupled cavity mode is visible.Fig. 5(c) shows RC spectra taken at ∼ 6.5°around the anticrossing at different laser fluences.The collapse of the two polariton peaks into one peak signifying the transition to the weak coupling regime is observed above 25 µJ cm −2 .The LPB and UPB energies extracted using the coupled oscillator model (Supplementary Figure S5) are shown in Fig. 5(d).As the polariton density is increased, the LPB and UPB approach each other almost symmetrically, converging to the exciton energy.The corresponding normalized Rabi splitting (Ω/Ω max , where Ω max is measured at low fluence) are shown in Fig. 5(d,e) as a function of the total polariton density.
In this experiment, the cavity mode is considerably above the X A energy, which therefore is not coupled to the cavity.Hence, the extracted Rabi splittings are fitted with a theoretically predicted trend of Ω for the NB excitation regime (Supplementary Note S8).A nonlinear polariton coefficient β = 0.86 µeVµm 2 is extracted by differentiating the fitted function with respect to the polariton density.Comparing our results to X A intralayer-exciton-polaritons in monolayers in similar cavities [16], we observe that the nonlinearity coefficient for dipolar interlayer polaritons is about an order of magnitude larger.This is in a good agreement with the theoretically predicted intrinsic nonlinearity of hybridized interlayer polaritons (Supplementary Note S8), and with our experimental data comparing hIX and monolayer X A outside the cavity (Supplementary Note S7).
In summary, we report the nonlinear exciton and exciton-polariton behaviour in MoS 2 homobilayers, a unique system where hybridized interlayer exciton states can be realized having a large oscillator strength.We find that nonlinearity in MoS 2 bilayers can be enhanced when both the intralayer and interlayer states are excited simultaneously, the regime that qualitatively changes the exciton-exciton interaction through the hole crowiding effect introduced theoretically in our work.In this broad-band excitation regime, the bleaching of the hIX absorption occurs at 8 times lower hIX densities compared to the case when the interlayer excitons are generated on their own.In addition to this, we find that the dipolar nature of hIX states in MoS 2 homobilayers already results in 10 times stronger nonlinearity compared with the intralayer excitons in MoS 2 monolayers.Thus, we report on an overall enhancement of the nonlinearity by nearly two orders of magnitude.Thanks to the large oscillator strength, hIX can enter the strong coupling regime in MoS 2 bilayers placed in microcavities, as realized in our work.Similarly to hIX states themselves, dipolar polaritons also show 10 times stronger nonlinearity compared with excitonpolaritons in MoS 2 monolayers.We expect that in microcavities where the cavity mode is coupled to both hIX and X A in MoS 2 bilayers, and the excitation similar to the broad-band regime can thus be realized, the nonlinear polariton coefficient will be dramatically enhanced owing to the hole crowding effect, allowing highly nonlinear polariton system to be realized.We thus predict that MoS 2 bilayers will be an attractive platform for realization of quantum-correlated polaritons with applications in polariton logic networks [20] and polariton blockade [21,22].

METHODS
The hBN/MoS 2 /hBN heterostructures were assembled using a PDMS polymer stamp method.The PMMA spacer for the monolithic cavity was deposited using a spin-coating technique, while a silver mirror of 45 nm was thermally evaporated on top of it.
Broad-band excitation was used to measure the reflectance contrast (RC) spectra of the devices at cryogenic temperatures (4K), defined as RC = (R sub − R BL )/R sub , where R sub and R BL are the substrate and MoS 2 bilayer reflectivity, respectively.
For the magnetic field studies the same RC measurements were performed using unpolarized light in excitation with polarizers, λ/4 polarizers and λ/2 waveplates in collection, to resolve σ + and σ − polarization.The low temperature measurements using the tunable cavity were carried out in a liquid helium bath cryostat (T=4.2K)equipped with a superconducting magnet and free beam optical access.We used a white light LED as a source.RC spectra were measured at each ∆L and are integrated over the angles within 5 degrees from normal incidence.The RC spectra measured in the cavity are fitted using Lorenzians.The peak positions are then used to fit to a coupled oscillator model, producing the Rabi splitting and the exciton and cavity mode energies.
The measurements on the monolithic cavity were performed in a closed loop helium flow cryostat (T=6K).For the power-dependent RC experiments, we used supercontinuum radiation produced by 100 fs Ti:Sapphire laser pulses at 2 kHz repetition rate at 1.55 eV propagating through a thin sapphire crystal.The supercontinuun radiation was then filtered to produce the desired narrow-band excitation.
All the exciton and polariton densities were calculated following the procedure introduced by L. Zhang et al. [16], taking into account the spectral overlap of the spectrum of the excitation laser and the investigated exciton peak (see further details in Supplementary Note S4).
In this section of Supplemental Materials we describe details of theoretical description of MoS 2 homobilayers.Properties of excitons in a bilayer system have been discussed in the main text, and we support them by modelling.The optical response of the system is characterised by the response of three different species of quasi-particles, namely X A , hIX and hX B (see main text).Here we provide an intuitive picture of the homobilayer physics and estimate the system parameters.
In particular, we estimate exciton Bohr radii and a hole tunnelling rate as relevant parameters when studying nonlinear properties.
We consider a MoS 2 homobilayer system with 2H stacking (see Fig. S1).This is comprised of two parallel layers of MoS 2 , with centres located at a distance d from each other (charge separation distance).The physics of bilayers is defined by properties of electrons and holes that interact through the Keldysh-Rytova potential Supplementary Figure S1.Side view of 2H-stacked MoS2 bilayer.Blue spheres are Mo (molybdenum) atoms and green are S (sulphide) atoms.We picture a hole wave function in each layer, constructing a hybrid state through tunneling process (hole delocalisation).This allows for hybrid hXB and hIX excitons.We define the distance between the layer as the distance between particle centres of charge.In absence of external fields, the centre of charge is located in between the S planes [44].[45][46][47], being different for in-plane and out-of-plane interaction [48].Within the k • p framework, electrons and holes are treated as particles with effective mass provided by a band dispersion.In MoS 2 the typical values for effective masses are 0.46 m e for conduction bands and 0.56 m e for the valence bands (with m e being the electron mass) [11,49].The attractive Kledysh-Rytova potential has a different form depending on the relative position between particles.We call V intra KR the attractive potential of particles being in the same layer, and V inter KR the attractive potential of particles being in separate layers.In momentum space the different potentials read as where e is the electron charge, 0 is the vacuum permittivity, in an average environment permittivity, r 0 is a screening length (defined as for monolayers), and q is an exchanged particle momentum [48].We compute a binding energy of an exciton bound state by assuming the Keldysh-Rytova attractive potential and free particle dispersion defined by the effective masses.To provide a simple understanding of the system, we approach the problem using an ansatz wavefunction φ(ρ) = 2/πα 2 exp(−ρ/α), with ρ being the in-plane projection of electron-hole distance and α the exciton Bohr radius.This describes well an internal structure of an exciton, and gives the information about its shape, collected in the Bohr radius.
In Fourier space, the function φ(q) reads and hX B , respectively.We considered a band-gap of 2.12 eV, and spin-orbit splitting of 13 meV for conduction band and 150 meV for valence band, as suggested by the ab initial calculations [49].Here we estimate the interlayer distance d thanks to experimental knowledge of energy distance between exciton energy modes.Note the peak of X A mode is not shifted by any tunneling, while X B and IX are coupled through tunneling constant J [50].That is, the theoretical energy distance ∆ 0 between the two modes has to be fixed to match the experimental result ∆ = 109 meV.With considering two coupled harmonic oscillators, we find the corrected energy shift to be By matching with the experimental data for the energy distance between modes, we Bohr radius grows with distance due to the reduced attraction between particles in separate layers.On the contrary, by increasing the interlayer distance, we see the reduced screening for particles in the same layer, resulting in a decrease of the Bohr radius and consequent increase of the binding energy.

SUPPLEMENTARY NOTE S2: COUPLED OSCILLATOR MODEL
A full picture of our system, corresponding to MoS 2 homobilayer, has to take into account 4 different modes coupling to each other: only IX and X B hybridise through a tunneling parameter, while X A and X B can couple with the cavity due to their high oscillator strength.We can then simplify this picture by rewriting the IX and X B states in terms of the new basis of hybridised modes, hIX and hX B , as defined in the main text, all of them now capable of a coupling with the cavity mode.The corresponding Hamiltonian reads where E c in the energy of the cavity mode, and E X A , E hIX , E hX B denote energies of the respective excitonic modes.Here, Ω X A , Ω hIX , and Ω hX B are corresponding matrix elements for light-matter coupling (Rabi splittings).
Due to the large energy separation between the resonances, each splitting can be fitted to a two level oscillator model independently.In our case, the spectra from the open cavity scans at piezo voltages close to resonant anticrossings between the cavity and either X A or hIX, were fitted with Lorenzian functions.The results were then fitted to the respective Hamiltonians of two coupled oscillators, such that  and the values of Rabi splittings and resonant energies were extracted.The results of the fit are shown in Fig. S3 In the presence of an out of plane magnetic field of 8 T, the same scans were repeated with unpolarised excitation and detecting opposite circularly polarized   of the Brillouin zone [52,53], while the slow component is possibly related to radiative [54] or defect-mediated non-radiative recombination [55].We conclude that in our MoS 2 bilayer sample the fast and slow decay times for both X A and hIX are much longer than the temporal width of the probe pulses (≈ 150 fs).
As there are many decay and dephasing the full treatment of possible processes if formidable.Here, we address as the key effect the decay process due to Coulomb scattering [40].As a result, the Coulomb scattering induced decay is proportional to the Coulomb scattering matrix where E is the energy of involved particles and V dir,exch (q) are direct and exchange particle scattering amplitudes, as discussed in Sec.S9.Note that in the main text for brevity we refer to the combined effect of different Coulomb-based processes using the combined interaction constant V Coul .The discussed decay properties of generic particles combine into the final shape of spectral peak L(E), see Fig. 4(a) [57]: where E 0 is the position of the peak, Γ R is the radiative decay rate and Γ NR is the non-radiative rate.Analysing the experimental data with minimal square method, we estimate the non-radiative decay rate to be of the order of ∼ 1 meV for direct excitons and ∼ 10 meV for indirect excitons.The ratio between the two agrees well with the estimates for the interaction constants (see the discussion below).With combining both radiative and non-radiative bleaching, we find the experimental data are well described by the relation where the discussed experimental observation are well described by ξ NR ∼ 10, and ξ sat ∼ 7 describes the nonlinear saturation of the Rabi splitting.The origin of the saturation term comes from the interlayer exciton phase space filling, and is reminiscent to phase space filling effects discussed in Ref. [29].
Finally, let us consider the exciton-exciton Coulomb scattering.We provide here the estimates interaction constants, with considering our gained knowledge of the effective interlayer distance d and particles Bohr radii α D,I .We build the interaction by following the procedure in Refs.[19,58].As a signature of particle indistinguishability, we have to consider both direct scattering V dir (q) and exchange scattering V exch (q) amplitudes, with q being the exchanged momentum.We omitted particle momenta as we consider total momentum to be zero.Explicitly, the two contributions read and where we define Φ tot (r e , r h , r e , r h ) as the sum of mutual particle interaction.For the scattering elements above, we shall consider two separate cases for direct and indirect excitons, as both wavefunctions and Coulomb terms differ.
With direct excitons, we have V dir X−X (q = 0) = 0, and We plot characteristic scattering amplitudes for direct excitons, where the direct process is zero at the negligible exchanged momentum (orange line), and the blue line corresponds to the exchange processes.
A is the sample area.Finally, the indirect exciton exchange potential has to be evaluated as Figs. S10(a) and (b) show the dependence of scattering amplitudes on the interlayer distance.We note the behavior of exchange scattering amplitude for indirect excitons, becoming negative past threshold distance (see Fig. S10(a)).Characteristically, for intralayer (i.e.direct) excitons we have a zero direct scattering amplitude, with non-zero contributions only due to the particle exchange process.
With the estimated parameters, we find interlayer exciton-exciton interaction to have a scattering constant of V I−I (q = 0) 2.5 µeV µm 2 , setting the scale for the Coulomb-based interactions.

SUPPLEMENTARY NOTE S9: THEORY FOR ENERGY SHIFT
The out-of-cavity results (see Fig. 3 main text) reveal an opposite behavior of the sample response depending on the excitation regimes, where narrow bandwidth (NB) and the broad bandwidth (BB) regimes are considered.In the NB case, we observe a blueshift for the hIX peak, as already reported in literature for dipolar excitons [19,59,60], while in the BB regime we observe a redshifted signal.The latter is unexpected, as in the system with dipolar excitons and predominantly positive scattering matrix elements for exchange terms, the emergence of some effective attractive nonlinearities requires a new possible mechanism.Below, we motivate the emergence of such a mechanism unique to the bilayer system.
First, let us analyse the energy of system as the expectation value of the full Hamiltonian Ĥ.We consider electrons and holes distributed over bilayer as shown in Fig. 1(d k,k ,q We define the exciton creation operators in the form By following a procedure for accounting non-bosonic correction at increasing order [43], we first commute the Hamiltonian with the product of exciton operators composite excitons.This reveals three possible creation potentials V µ,ν D−D , V µ,ν I−I , and V µ,ν D−I , which noticeably generate cross-flavour interaction.The presence of creation potentials, as well as quadratic scaling of these terms in the total energy, is the signature of nonlinear behaviour.Terms in lines 2 and 3 of Eq. (S22) are the energy shifts for direct and indirect excitons due to phase space filling within the same exciton flavour.The last term in line 4 is the cross interaction of direct and indirect excitons, allowing for extra energy shifts in modes that are not statistically independent.Namely, the corresponding operators for intralayer and interlayer excitons do not commute as they share a hole, but formed by different electrons in conduction bands.This leads to the negative valued commutator, and here we identify the origin of unusual redshift, as observed in the work.We refer to this effect as hole crowding, where hole population being shared between the two different particles (hIX and X A ).This is different from statistical deviation of same excitons, where both carriers are exchanged, which we simply refer as interlayer exciton phase space filling, in analogy to the monolayer case.The effect of intra-flavour terms has already been discussed in Ref. [29], and leads to positive energy shifts.However, the cross-flavour terms only emerge in bilayers, and to date remained unexplored.
We evaluate the considered term in q 0, as it gives the dominant contribution

FIG. 1 .
FIG. 1. Homobilayer MoS 2 and its optical response.a, Bright field microscope image of an encapsulated BL MoS 2 transferred on top of a DBR.Scale bar: 10 µm.b, Schematic side-view of the fabricated heterostructure comprising a BL MoS 2 sandwiched between few-layer hBN.c, Reflectance contrast (RC) spectrum of the sample measured at low temperature (4 K) showing three distinct

FIG. 2 .
FIG. 2. Strong exciton-photon coupling in MoS 2 bilayers.a, Schematics of the tunable open microcavity composed of a bottom DBR and a top semi-transparent silver mirror.b, c, Low temperature (4K) RC spectra measured as a function of the cavity-exciton detuning (∆ = E cav −E exc ) for cavity scans across b X A and c hIX energies.White dotted lines show the fitting obtained using the coupled-oscillator model providing the Rabi splittings Ω hiX = 19 meV and Ω X A = 38 meV.d, RC spectra measured for the cavity-exciton detunings in the vicinity of the anticrossing between hIX and the cavity mode.e, Dipolariton dispersion measured with circularly polarized detection for 8 T magnetic field.The orange and black solid curves are the coupled oscillator model fits for σ + and σ − detection, respectively.The positions of the Zeeman-split hIX peaks are shown by dashed lines.f, σ + (orange) and σ − (black) RC spectra measured at 8 T at the hIX-cavity anticrossing.Fitting with two Lorentzians (solid lines) is shown.

FIG. 3 .
FIG. 3. Exciton nonlinearity in MoS 2 bilayers.a, b, c, RC spectra measured with the NB (FWHM=28nm) excitation for the X A (a) and hIX (b), and with the BB (FWHM=50nm) excitation (c) at different fluences.The dashed curves are guide for the eye.d, e, The energy shift ∆E (top), linewidth variation ∆FWHM (middle) and normalized integrated intensity (bottom) as a function of the exciton density for the hIX (d) and X A (e). Solid (open) symbols show the results for the BB (NB) excitation.For the normalized intensity we divide the intergrated intensity at each laser power by that at the maximum intensity.

FIG. 4 .
FIG. 4. Theoretical model for nonlinear optical response in MoS 2 bilayers.a, b, Schematic diagram showing exciton generation under the NB (a) and BB (b) excitation.In (a) only generation of hIX is shown.In (b), the holes of the two excitonic species share the same valence band.c, Theoretically calculated absorption spectra for the BB excitation case (see Supplementary Note S8), providing qualitative agreement with the experiment.The dashed black curves are guides for the eye.

X
A .Using the estimated nonlinear coefficients caused by the hole crowding, we model the RC in the BB regime and qualitatively reproduce the strong bleaching and redshift for hIX at the increased density.

FIG. 5 .
FIG. 5. Nonlinear behaviour of dipolaritons.a, b, Reflectance contrast spectra measured at different laser fluences for the MoS2 bilayer placed in a monolithic cavity.(a) The low fluence case (0.6 µJ cm −2 ).A clear anticrossing at 6.5°is observed.Dashed red lines show the results of the fitting using a coupled oscillator model, with two polariton branches LPB and UPB formed.White and orange lines show the energies of the uncoupled cavity mode and hIX state, respectively.The vertical line marks the anticrossing angle.(b) The high fluence case (58.5 µJ cm −2 ).A complete collapse of the strong coupling regime is observed, with the disappearance of the anticrossing and transition into the weak coupling regime.c, RC spectra measured at the anticrossing at 6.5°as a function of the laser fluence.d, Measured UPB and LPB peak energies at 6.5°as a function of the laser fluence (see top axis) and the corresponding polariton density (bottom axis).e, Symbols show the Rabi splittings normalized by the Rabi splitting measured at the lowest power (Ω/Ωmax) as deduced from d,.The line shows the fitting using our theoretical model (Supplementary Note S8).

)
With the given wave function [Eq.(S2)], and the given potentials [Eq.(S1)], we can find the binding energy and the Bohr radius of the excitons by minimizing the energy of the system.This procedure is performed in the range of possible interlayer separation d.Fig.S2(a)shows the energy change with the separation.The blue curve corresponds to the hIX mode, while orange and red curves correspond to X A

Supplementary
Figure S2.Evolution of particle properties with interlayer distance.a) Energy of quasiparticle modes with distance.The blue line is the hIX mode, and the orange and red respectively are XA and hXB.The redshift of direct modes with the increase of interlayer distance is due to reduced screening effects, as characteristic of Keldysh-Rytova potential in Eq. (S1).On the contrary, the hIX mode sees a blueshift that is due to the lower attraction of particles located in separate layers.The purple dashed line is the distance to match with the experimental data.b) Dependence of Bohr radius with interlayer distance.Here we only report one direct exciton as both A and B excitons have the same behavior.The orange one is the evolution for direct excitons, and the blue line is for indirect excitons.Both plots (a) and (b) show similiar increase (or increase) with distance as the screening mechanism affects the considered parameters the same way.This analysis reveals the hIX Bohr radius to be roughly twice as much as the X one.extract d = 5.5 Å and estimate the effective tunneling rate to be J = 45 meV.As a consequence, 21% of X B oscillator strength is transferred to IX mode, in agreement with previous observations [11].Finally, Fig. S2(b) shows the evolution of particle Bohr radii with the interlayer distance.We respectively call the Bohr radius of direct and indirect exciton α D and α I .The blue curve is the Bohr radius of hIX, while the orange one described the X modes.With X A and hX B being very similar, we describe both with one orange curve.The energy separation is provided only by the spin-orbit splitting.Typical values of α D are approximately 1 nm, with α I being approximately 2 nm.Note the opposite behavior of α D and α I with distance.hIX

Supplementary
Figure S5.Coupled oscillator model fits on monolithic cavity.Monolithic cavity dispersion near hIX with couple oscillator model fits.Blue and green points show the extracted peak positions of the LPB and UPB, respectively, from the spectra at each angle.These are then fitted to the two level coupled oscillator model and the solutions LPB and UPB are shown as dark orange dashed curves.The extracted cavity mode and exciton energy are shown as white and light orange dashed curves, respectively.consist of a white-light continuum (WLC), generated from a 1-mm thick sapphire plate pumped by focusing the 800 nm output of the main laser.Pump and probe pulses are synchronized by means of a motorized delay stage.Pump and probe beams are then collinearly combined by a thin dichroic beam splitter and focused on the sample using an objective lens (NA=0.3),resulting in a 4µm spot size.The samples are placed in a closed-cycle helium cryostat reaching a temperature of 6 K.The spatial overlap of the sample with the pump and probe spots is obtained by a three-axis (xyz) mechanical translation stage coupled to a home-built imaging system.After the interaction with the excited sample, the reflected probe pulse is collimated by the same objective lens and then sent to a spectrometer equipped with an electronically cooled Si CCD, to measure the differential reflection (∆R/R) Supplementary Figure S6.a) Schematics of the pump-probe microscopy setup used for the experiments on the MoS2 bilayers.b, c) Transient reflectivity traces for XA (b) and hIX (c) taken at 643 nm and 622 nm respectively.The red curves refer to the fitted bi-exponential decay function.

Figure
FigureS6 (b,c) shows the temporal exciton dynamics measured as the absorption bleaching signal in transient reflectivity, taken at the peak wavelengths of the pumpprobe traces, 643 nm and 622 nm for X A and hIX respectively.For this experiment ) where δ(x, y, θ) = x 2 + y 2 − 2xy cos θ.With indirect excitons, the direct scattering amplitude recovers the capacitor formula V dir I−I (q = 0) = e 2 /( 0 A)d, where Supplementary Figure S10.Scattering amplitudes for exciton-exciton Coulomb interaction.(a,b) We show the change of scattering amplitudes with the interlayer distance.In (a), note the behavior of exchange scattering amplitude for indirect excitons, becoming negative past a threshold distance.Direct scattering amplitude grows linearly with distance.(b) ) [main text].We call N I(D) the number of indirect (direct) excitons in the sample, and n I(D) the particle density.The spin indices are omitted for brevity.The Hamiltonian of the system can be written asĤ = ĤT + Ĥsm + Ĥdf , (S14)where ĤT is the kinetic term, Ĥsm and Ĥdf are the Coulomb interactions.With Ĥsm we refer to interacting particles belonging to a same band, and Ĥdf corresponds to different dispersion bands.Explicitly, the kinetic energy readsĤT = k ε t c (k)â † k âk + ε t v (k)b † k b k + ε b c (k)c † k c k ,(S15)where â † k and b † k are creation operators for conduction and valence bands of the top layer, respectively.ĉ † k is an electron annihilation operator for the conduction band of the bottom layer.Each operator is labelled with a crystal momentum k. ε t(b) c(v) (k) are dispersions for conduction (valence) bands in the top (bottom) layer.The interactions are mediated through the Keldysh-Rytova potential.The corresponding â † k+γ e P ĉk−γ h P , (S19) with D † µ (Q) [ Î † ν (P)] and φ µ D (k) [φ ν I (k)] being the direct [indirect] exciton creation operator and wave function, respectfully.Q, P are crystal momenta, and µ, ν are state indices.When we omit the µ and ν indices, and the total momentum, we refer the ground state at crystal momentum Q = 0. To take into account for particle non-bosonicity and consequent nonlinear behaviour, we consider the expectation value of a system over a multi-particle state created by exciting the vacuum state |Ω 0 .It includes N D direct excitons and N I indirect excitons.The expectation Ω 0 | DN D ÎN I Ĥ D †N D Î †N I |Ω 0 then denotes the total energy of the many-body system.
observed in experiments, and we consider particles being in the ground state as a main occupation at low temperatures.With this, we rewrite the mutual energyshift ∆E I = µ,ν,q V µ,ν D−I (q) Ω 0 | DN D ÎN I D †N D −1 D † µ (q) Î † ν (−q) Î †N I −1 |Ω 0 as ∆E I = −ξ e 2 4π 0 A Iα, (S23)with A being the sample area, and parameter ξ of the order of unity, ξ ∼ 1, is tuned to match with experimental data setting an effective area.Here, α = (α −1 D + α −1 I ) −1