{'fraction_non_alphanumeric': 0.028864012850760505, 'fraction_numerical': 0.011696199989960344, 'mean_word_length': 4.533888888888889, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Highly non-linear interlayer exciton-polaritons in bilayer MoS2

Biswajit Datta, Mandeep Khatoniar, 2 Prathmesh Deshmukh, 2 Rezlind Bushati, 2 Simone De Liberato, Stephane Kena Cohen, and Vinod M. Menon 2, a) Department of Physics, City College of New York, New York, NY, USA Department of Physics, Graduate Center of the City University of New York (CUNY), New York, NY, USA School of Physics and Astronomy, University of Southampton, Southampton, UK Department of Engineering Physics, École Polytechnique de Montréal, Montréal, Quebec, Canada

Realizing nonlinear optical response in the low photon density limit in solid state systems has been a long-standing challenge. Semiconductor microcavities in the strong coupling regime hosting light-matter quasiparticles called exciton-polaritons have emerged as an attractive candidate in this context. However, the weak interaction between these quasiparticles has been a hurdle in this quest. Two-dimensional transition Photons are becoming ubiquitous in emerging quantum technologies like quantum communication and metrology due to the ability to propagate long distances while being robust against decoherence 1,2 . To further extend the span of their utilities, platforms to achieve and implement strongly interacting photons in solid-state systems are much-sought after 3 .
Conventional materials do not exhibit nonlinear response at power levels associated with single photons. In this context, remarkable advances have been made in cold atom systems to realize interacting photons at the single-particle level 4 . A strong contender for the generation of strongly interacting photons in the solid-state are exciton-polaritons formed via non-perturbative coupling of cavity photons with excitonic resonances. Although they can be modeled as non-interacting quantum fluids at low densities, beyond a critical density, saturation and short-range exchange interactions become significant and give rise to various phenomena like the appearance of spontaneous coherence, parametric down-conversion, and superfluidity 5 . Strong spatial confinement along with Coulomb interactions give rise to even stronger correlations between polaritons that can no longer be described using a mean-field theory. Such interactions give rise to non-Poissonian statistics of laser transmission, dubbed as polariton blockade. Preliminary evidence of such non-classical correlation was recently observed in confined polariton systems in GaAs 6,7 . The small ratio of interaction to dissipation rates resulted in only a weak violation of classical correlations. Formation of polaritons with excitons that possess a permanent dipole moment has been shown to enhance the polariton interactions both in resonant 8,9 and non-resonant excitation schemes 10 .
Transition metal dichalcogenides (TMDCs) have garnered much attention for their exceptional optoelectronic properties and have been used to demonstrate wide array of fundamental phenomena as well as technological applications 11 . Furthermore, their two-dimensional (2D) nature allows for the formation of heterojunctions or homojunctions with arbitrary twist angles resulting in emergent properties [12][13][14][15] . Their strong binding energy along with large oscillator strength, renders them capable of forming polaritons at room temperature 16 , which has also been shown to retain the intriguing properties of the 2D excitons [17][18][19] . Polariton interactions in TMDCs under different configurations are being studied extensively.
Realization of Fermi polarons 20 , trion polaritons 21,22 , and excited-state exciton-polaritons 23 have shown great potential in their abilities to harness strong polariton interactions. Moiré exciton-polaritons were recently demonstrated using heterobilayers of WS 2 /MoSe 2 where the electronic confinement potential arising from the twisted heterostructure was shown to enhance the nonlinearity 12 . However, all the platforms above rely on short-range exchange interactions that pose a bottleneck in realizing few polariton non-linearity under current experimental capabilities.
In this context, a very attractive possibility is the use of spatially indirect interlayer excitons (IE) in TMDC heterostructures 24 that can have a permanent dipole moment and hence support highly interacting dipolar polaritons. However, their inherently low oscillator strength creates an impediment in reaching the strong coupling regime. Furthermore, the large in-built interfacial electric fields make electrical tuning of the energies of the IE far more difficult. Bilayer MoS 2 provides a highly attractive platform to solve these issues [25][26][27][28] .
The IE in naturally stacked 2H bilayer MoS 2 remarkably has both an out of plane (static) dipole moment and an in-plane (oscillating) dipole moment. Due to the in-plane dipole moment, the IE in MoS 2 bilayer has an oscillator strength of approximately 36% of that of the intralayer A exciton along with strong absorption that is visible up to room temperatures.
Moreover, a strong response to DC electric fields has been demonstrated in these systems, thus confirming their dipolar nature 25,26,28 . In this work, we achieve strong coupling of microcavity photons with the IEs (along with intralayer A and B excitons) in bilayer MoS 2 .
The IE polariton shows 10 fold enhancement of the polariton nonlinearity compared to the intralayer A exciton, making it an appealing platform to realize strongly interacting polaritons in condensed matter systems. In addition, the ease of fabrication and realization of multi-polariton species in this system makes it a practical and fundamentally interesting material for studying polariton physics. IE 1 and IE 2 excitons in which electrons are in layer 1 and layer 2 respectively remain degenerate at zero external bias. The hole is always delocalized among both layers. The black arrows connect the bands that form the IE. This peculiar charge distribution is at the heart of producing both in-plane and out-of-plane dipole moments of IE. Fig. 1b shows the charge distribution in real space for IE 1 and IE 2 . Owing to the spatial separation of the electron and the hole, it acquires a permanent dipole moment, the nature of which has been a topic of much interest recently [25][26][27][28] . IEs in MoS 2 homobilayer can be thought of as an admixture of B intralayer exciton with an optically dark but electric field tunable charge-transfer exciton, which is typically found in TMDC heterobilayers 29 . As a result, IEs in MoS 2 homobilayer acquires both a strong oscillator strength and electric field tunability 29 . Fig. 1c shows the white light differential reflection of the bilayer sample at 77 K. The peaks at 1.933 eV, 1.995 eV, and 2.10 eV correspond to the A, IE, and B excitons, respectively. The absorption of the IE is found to be 36% of the A exciton's absorption. Interestingly, the absorption of the IE is discernible even at room temperature, which provides a straightforward method to identify the MoS 2 bilayer flakes after mechanical exfoliation. Supplementary Note II for details on the coupled oscillator model. We notice that the Rabi splitting of the IE polariton remains observable even at 77 K but not at room temperature due to increased damping at a higher temperature. On the other hand, 1s state of A and B excitons remain strongly coupled even at room temperature due to their larger oscillator strength. See Supplementary Note III for the data at 77 K.
We measured the white light reflectivity around the IE and A resonance and monitor the energy shifts of the polariton branches, see Supplementary Note V for the details. Fig. 3a shows the energy of the upper and lower branch of the IE polariton as a function of polariton density at the inplane wavevector corresponding to the zero detuning k (where the cavity mode and exciton energy are degenerate). We fit Lorentzians at each density to obtain the energy of the upper and lower polaritons. See Supplementary Note IV for details of the polariton density calculation and Supplementary Note V for the raw data of the density dependent differential reflection data of IE polariton. We can see from Fig. 3a that the lower branch of the IE polariton moves with power more than the upper branch at the zero-detuning k || . We also observe a simultaneous increase of the zero-detuning k || and the decrease of Rabi splitting (at the zero-detuning k || ) with increasing density. This suggests the presence of both exciton energy renormalization and saturation effect. See Supplementary Note VIII for the density dependent Rabi splitting of IE polariton. Taking the derivative of the E LP with respect to density, we calculate the strength of the nonlinearity, g LP . Since the rate of blue shift saturates at high density (olive color curve in Fig. 3a), g IE LP reduces with polariton density as seen in Fig. 3b. Note that at the lowest powers accessible in our measurement, the IE LP already shows a power-dependent blueshift. On the other hand, within the accessible range of the polariton density in our experiment, g A LP remains nearly a constant (see Supplementary Note IX and Supplementary Fig.S8 for the blueshift data of the lower branch of A polariton). The obtained strength of the non-linearity in the low density limit for IE and A polaritons are g IE LP ∼ (100 ± 2) µeV µm 2 and g A LP ∼ (10 ± 0.2) µeV µm 2 respectively. For IE polariton, we find a significant ∼10 fold increase in the nonlinearity compared to the A exciton-polariton. This g A LP value for pol-1 is comparable to the previous report for A exciton-polariton in monolayer WSe 2 23 .
We now proceed to discern the contributions from nonlinearities arising from renormalization of the exciton energy (∆E exc−exc ) and reduction of the oscillator strength of the exciton due to phase space filling (∆E sat ) to the energy shift of pol-2 and pol-3. Fig. 4a shows the Hopfield coefficients of all the excitons and the cavity as a function of k which is later used to plot ∆E IE LP as a function of Hopfield cavity fraction (C). The Hopfield coefficients of A exciton-polariton is given in Supplementary Note VII. We measure the energy shift ∆E IE LP (of pol-2) as a function of the Hopfield cavity fraction, see Fig. 4b. The interaction strength due to phase space filling in a single excitonic system can be written in terms of g LP = 4g SAT |C||X| 3 , where X is the Hopfield coefficient for the exciton. Although this is a multi exciton system, the contribution of the other excitons close to the zero detuning k of the IE is small. Neglecting other excitons we can write g IE LP = 4g SAT |C|(1 − |C| 2 ) 3/2 which is a non-monotonic function of the cavity Hopfield coefficient as shown in the Fig. 4b inset (red curve). Similarly, considering only exciton-exciton interaction, the energy scales as g LP = g XX |X| 4 = g XX (1 − |C| 2 ) 2 which monotonically decreases with increasing |C| (cyan curve in Fig. 4b inset). We find that the measured nonlinear response of the IE polariton (at the Rabi splitting of pol-2 and pol-3) is the highest at an in-plane wave vector larger than zero detuning k , as seen in Fig. 4b. We also notice that ∆E IE LP changes non-monotonically with the cavity Hopfield coefficient. The qualitative resemblance of the measured ∆E IE LP with the non-monotonic saturation effect suggests that it is the dominant mechanism dictating the nonlinearity in our system. However, the appearance of the ∆E IE LP peak at a |k | > zero detuning k and increase of zero detuning k at higher density points to a sizable contribution of the ∆E exc−exc . Fig. 4c shows a schematic of the polariton energies at two different powers for ∆E sat > ∆E exc−exc . Due to both the blue shift of the exciton and saturation effect, the lower branch moves with increasing power more than the upper branch. Analyzing the blueshift of the IE lower polariton branch and redshift of the IE upper polariton branch we calculate the magnitude of ∆E exc−exc and ∆E sat in our system. At the zero detuning k we write the energy shift of pol-2 as ∆E IE LP = ∆E exc−exc + ∆E sat and the energy shift of pol-3 as ∆E IE UP = ∆E exc−exc − ∆E sat . From these equations we obtain ∆E exc−exc and ∆E sat as a function of polariton density, as shown in Fig. 4d. This shows magnitude of ∆E sat in our system is larger than ∆E exc−exc by roughly a factor of 2.
The polaritons interact with each other due to their excitonic part in the wavefunction.
In general, the Coulomb interaction of excitons has both direct and exchange contributions.
For a regular intralayer exciton (such as A exciton in our system), the direct Coulomb energy goes to zero for small wavevector 30 . Hence for intralayer excitons in TMDC, the exchange term dominates. On the contrary, due to the separation of electrons and holes in two different layers for interlayer excitons in heterobilayer TMDC, the exchange term is negligible. At the same time, the direct term is non-zero even at small wavevectors for interlayer excitons 30 . As a result, the direct Coulomb term is the dominating energy correction for interlayer excitons as shown theoretically in ref. 30 . The direct interaction for small momenta scales with the fourth power of excitonic Bohr radius 30 . Moreover, the saturation density is inversely proportional to its second power and depends on the details of the excitonic wavefunction 31 . A smaller binding energy has been reported in the literature for IE compared to the A exciton (112 meV and 163 meV respectively 26 ), as expected from the increased real-space charge separation in the former. The corresponding larger Bohr radius will thus lead to higher interaction energy and smaller saturation density. This is a plausible reason for our observed 10 fold enhancement in the nonlinearity for the IE compared   µeV µm 2 is almost 10 times larger than the g LP of A exciton polariton ∼ (10 ± 0.2) µeV µm 2 .
The density error bars take into account the error in measuring power and determining the energy of the polariton and linewidth. The error bars in g LP take into account the error in determining density and the energy of the polariton.   Angle is converted to cavity Hopfield coefficient using the previous plot. We note that the blueshift of the lower branch of IE (∆E IE LP ) is highest close to the zero detuning k , but the shape of the curve looks like a tilted parabola. This can be explained by taking both exciton-exciton interaction and saturation effect into account. Inset shows exciton-exciton interaction monotonically decreases with increasing |C| (cyan curve); on the other hand, the saturation effect initially rises with increasing |C| but then it decreases (red curve). It can be easily seen that superimposing these two graphs can explain the tilted parabola shape of ∆E IE LP in our data. Note that |C| is not 0.5 at zero detuning k because we do multiple exciton coupled oscillator fit. (c) Cartoon showing the IE polariton dispersion at two different excitation powers where the combined effect of exciton-exciton interaction and saturation effect is present. Here, X 1 and X 2 are the exciton energies at power P 1 and P 2 , respectively. In our system, even though the saturation effect outweighs exciton-exciton interaction, the zero detuning k increases with increasing power since the exciton blueshifts. We

II. DETAILS OF COUPLED OSCILLATOR MODEL
The five polariton branches in our experiment can be described with the eigen modes of the following five coupled oscillator model. 

IV. DETAILS OF THE POLARITON DENSITY CALCULATION
Polariton density is calculated from the Gross-Pitaevskii (GP) equation. The GP equation at wavevector k can be written as Here LP (k) = ω LP (k) is the energy of the lower polariton of an exciton and γ LP is the line width of the polariton. F P (k, t) is the coherent driving. The pulse width of our supercontinuum laser is 20 ps which is much larger than the polariton lifetime (30-50 fs).
This allows the pumping to be treated as continuous wave which is resonant with our cavity.
We can express the pumping term as F P (k, t) = F P (k)e −iω P t and ψ LP (k, t) =ψ LP (k)e −iω P t where ω P is the pump frequency. This simplifies the above equation to The F P (k) term can be written from the input-output relation.
Here C(k) is the photon Hopfield coefficient of the polariton branch. P int (k) is the incident power on the top mirror and η is the coupling coefficient. t top mirror is the transmission of top mirror, and τ trip is the photon trip time in the cavity. This gives η = 4.1 × 10 12 s −1 .
Equation S2 can be written as Polariton density at the wavevector k and frequency ω P is Since the polaritons have a finite width in energy, the total polariton density at the wavevector k can be found out by integrating over the all frequency range Here, ξ(k, ω P ) denotes the incident power density (power/frequency/wavevector).
The integral above can be numerically calculated but we can make an approximation to get a close form solution for the density. We note that the Lorentzian part of the integrand for γ LP → 0 is a Dirac delta function Using the above formula and taking the power density ξ(k, ω P ) constant over the narrow range of excitation we get Here, γ LP (k) is determined from the experimental data.
Total real space polariton density can be found out by summing ψ LP (k) 2 in the k-space for the experimental wavevector range If P 0 is the measured real space peak power density (power/energy) of the pulsed supercontinuum laser excitation then ξ = P 0 /N , where N is the number of k point in between k min and k max .   The upper polariton branch red shifts and lower polariton branch blue shifts with increasing density reducing the Rabi splitting. We also notice that the two branches move asymmetrically -the lower branch moves more than the upper branch at zero detuning k . This is because both exciton-exciton interaction and saturation give rise to blue shift for the lower branch but they produce opposite shift for the upper branch -as described in the main text.

IN OPEN CAVITY GEOMETRY
To examine the nonlinearity of the bare excitons even before they form polariton in a cavity, we studied power dependence on the same bilayer MoS 2 sample before depositing the top mirror. We observe that the amplitude of the differential reflection dip for both A and interlayer exciton reduces with increasing power; the FWHM of the lineshape also increases for both A and interlayer exciton with increasing power, see Fig. 4a and Fig. 4b. We analyze the area under the differential reflection curve as a function of excitation power. The area under the curve can be mapped to the oscillator strength assuming that the area under the differential reflection curve is proportional to the oscillator strength, see Fig. 4c and Fig. 4d. Here, f X and f X0 denote the oscillator strength of the exciton X at a given power value and at the starting power value respectively. We can see that the interlayer exciton absorption starts saturating at much lower power values than the A exciton. Comparing the power axis of Fig. 4c and Fig. 4d, we notice that the interlayer exciton reaches the same level of saturation at almost an order of magnitude lower power value than A exciton.   We plot the Rabi splitting of the IE polariton at the zero detuning k as a function of density, see Fig. S7. We mote that the usual formula for free 2D excitons Ω = Ω 0 √ 1+ n n sat does not fit our data well. The inadequacy of the free 2D exciton formula to describe the density dependence of the Rabi splitting was also observed in a recent study S12 on the Moiré exciton-polariton in a hetrobilayer TMDC.  S8. Energy of the lower branch of A exciton polariton as a function of polariton density at 7 K from which the g A LP is calculated and plotted in Fig.3b of the main manuscript.