Highly nonlinear trion-polaritons in a monolayer semiconductor

Highly nonlinear optical materials with strong effective photon-photon interactions are required for ultrafast and quantum optical signal processing circuitry. Here we report strong Kerr-like nonlinearities by employing efficient optical transitions of charged excitons (trions) observed in semiconducting transition metal dichalcogenides (TMDCs). By hybridising trions in monolayer MoSe2 at low electron densities with a microcavity mode, we realise trion-polaritons exhibiting significant energy shifts at small photon fluxes due to phase space filling. We find the ratio of trion- to neutral exciton–polariton interaction strength is in the range from 10 to 100 in TMDC materials and that trion-polariton nonlinearity is comparable to that in other polariton systems. The results are in good agreement with a theory accounting for the composite nature of excitons and trions and deviation of their statistics from that of ideal bosons and fermions. Our findings open a way to scalable quantum optics applications with TMDCs.

Supplementary Note 1 Experimental details

Parameters characterising strong coupling between photons and exciton and trion
To describe the polariton system where the photonic mode is strongly coupled to both excitons and trions we employ a model of three coupled oscillators, which can be described in a compact way by the matrix equation   E C (L) 1 2 Ω X where E C (L), E X , E T are the cavity photon, exciton, and trion energies; Ω X and Ω T are the cavity mode-exciton and the cavity mode-trion coupling strength. ψ is a three-component basis vector. In our system we can tune E C by changing the cavity length L. Eigenvalues of the matrix in the left-hand side of Supplementary Equation (1) give the polariton mode resonances of upper, middle, and lower polariton branches (UPB, MPB, and LPB, respectively). Their eigenvectors provide the Hopfield coefficients. We start by fitting this model to the experimental data in Fig. 1b (main text), and obtained the following values for the ground LG 00 mode: E X = 1646.0 ± 0.5 meV, E T = 1621.2 ± 0.5 meV, Ω X = 17.2 ± 0.5 meV, and Ω T = 5.8 ± 0.5 meV. Using these parameters we can plot polariton resonances and Hopfield coefficients for different values of cavity resonance energies E C (cavity lengths, L) expressed as a cavity-exciton detuning, δ C−X , in Supplementary Figure 1. We summarise Hopfield coefficients obtained for the four experimental cavity-exciton detunings used here in Supplementary Table I. Note, that these are the values for the low-density case. In the high density regime, as we discuss in the main text and the theory section, the nonlinearity quenches photon-trion coupling, and the polariton system can be described by the two-oscillator model.

Deduction of the total reservoir density in the case of nonlinear trion-polaritons
The scattering of polaritons with the excitonic disorder potential (which we estimate in the order of 10-15 meV from the inhomogeneous linewidths of exciton/trion emission) effectively populates the reservoir states 1-4 with a scattering rate of the order of polariton linewidth γ pol ∼ 3 − 5 meV. The polariton nonlinearity is driven by the total density of excitons and trions n tot excited in the system with a single pulse, i.e excitons and trions due to polariton density and due to reservoir. By measuring the total number of photons transmitted through the microcavity (MC) in a single pulse n phot (normalised to the cavity mode area) we can deduce this density as Here |X| 2 , |T | 2 , |C| 2 are exciton, trion and photon fractions of polariton. The bare cavity mode linewidth γ c ≈ 0.4 meV corresponds to the photon lifetime of τ c ≈ 4 ps. The third term in the above formula is dominant. It describes the ratio of absorbed to radiatively escaped polaritons as follows from a rate equation model, assuming that the lifetime of the reservoir is much longer than cavity photon lifetime and that there is no backscattering from the reservoir to polariton states. Indeed, suppose the density of polaritons excited with a short 100 fs pulse in the cavity is equal to N . This density may decay in the system by absorption due to scattering into the reservoir or by escape through the Bragg mirrors. The absorption rate into reservoir is given by dn r /dt = N/t r , where the scattering rate 1/t r = γ pol /h since γ pol γ c and n r is the density of polaritons absorbed into the reservoir. The escape rate into free space through  Bragg mirrors is given by dn phot /dt = N/t c , where 1/t c = γ c |C| 2 /h, and C is the photon fraction and n phot is the detected number of photons normalised to the cavity mode area A. Dividing the two equations above we obtain dn r /dn phot = t c /t r = γ pol /(γ c |C| 2 ) or dn r = γ pol γc|C| 2 dn phot . Integrating this equation over the duration of a single pulse we obtain the population of the reservoir with respect to the total number of the emitted photons within a single pulse n r = γ pol γc|C| 2 n phot , which corresponds to the third term in Supplementary Equation (2). There is a possibility that we may slightly overestimate the total excited exciton/trion density and hence underestimate polariton nonlineairity if there are backscattering processes from reservoir to polariton mode. On the other hand, since the absorption rate of the resonantly excited polaritons cannot occur at a rate larger than that given by the measured polariton linewidth γ pol we exclude an experimental underestimation of n tot and hence overestimation of polariton nonlinearity. For the case of resonantly driven trion-polaritons the reservoir and hence n tot consists of mostly trions, since the peak of exciton density of states is blue detuned by several ten's of meV from the trion.
Finally, we note some uncertainty of ∼ 20 % in the deduction of the total excited exciton/trion density may arise from the deduction of polariton linewidths of the transmission spectra. While trion-polariton spectra in Fig. 3 of the main text can be well fitted with Gaussian profile, the fitting of neutral exciton-polariton spectra Fig. 4a [main text] is not perfect, where the tails of the polariton spectra cannot be accounted for by a simple Gaussian form. This uncertainty in γ pol could lead to the ∼ 20 % uncertainty of the absolute values of the polariton nonlinear coefficients.
To deduce n tot , one has to measure the number of photons n phot emitted by the cavity in a single pulse after the resonant excitation with the 100-fs pulsed laser. To do that we shifted the position of the photon mode to the very negative detuning about 10 nm below the trion level and ramped up the power of the laser so that the power of the light transmitted through the microcavity is about 50 nW. The corresponding photon counts of the transmitted light on the CCD were measured as a reference value. By relating the measured photon counts of light emitted by polariton microcavity to this reference value it is possible to deduce the average power of the polariton emission in W, I pol . The density of photons emitted by microcavity in a single pulse is then calculated as n phot = I pol f ωA , where ω is the polariton energy, f = 10 3 Hz is the repetition rate of our pulsed laser and A ≈ 3 µm 2 is the mode area. Importantly, the very low repetition rate of our laser prevents the effect of dark exciton/trion reservoir excited in previous pulses.
The Hopfield coefficients used in Supplementary Equation (2) for the total reservoir density also have to be determined for each pump power, since the change of either cavity-trion or cavity-exciton Rabi-splitting would change these values. To illustrate this we fix the cavity-exciton Rabi-splitting and plot theoretical peak positions of LPB and MPB modes as a function of cavity-trion Rabi-splitting, which is shown in Supplementary Figure 2 for the cavity mode-exciton detuning used in Fig. 3 [main text]. The theoretical dependence of the photonic fraction of the polariton branches on the cavity-trion Rabi-splitting is shown in Supplementary Figure 3. Combining these two sets of data one can produce the "calibration" dependencies of the Hopfield coefficients. The photonic fractions |C| 2 vs MPB and LPB energies are shown in Supplementary Figure 4. This can be used to infer the values of the actual Hopfield coefficients of MPB branch in Fig. 3 of the main text for each pump power.
This method is valid since in the experiment cavity-trion Rabi-splitting is completely quenched before any effect is seen on the cavity-exciton Rabi-splitting, as discussed in the main text. Thus, to deduce total reservoir density we, first, determined the peak position of MPB mode by fitting each measured spectra, and then using this value we determined effective photon, exciton, and trion fractions of the MPB mode from the calibration curves. Finally, the total reservoir density for each pump power was calculated using Supplementary Equation (2).
3. Deduction of the total reservoir density in the case of nonlinear neutral exciton-polaritons Since the neutral exciton-polariton nonlinearities are weak in comparison to trion-polaritons they are studied at high pump when the strong coupling between photons and trions is quenched. In this case the polariton system can now be described by a simpler two-oscillator model, which results in two polariton branches UPB and LPB. ψ is a two-component basis vector. Separate MPB and LPB do not longer exist, since they recombine into a single MPB. This model is applicable for all polariton densities > 200 µm −2 , which is the case for all data shown in Fig Table II shows the values of the Hopfield coefficients for the same cavity mode-exciton detunings as in Supplementary Table I but now assuming that there is no trion-cavity mode coupling. Using the same approach as discussed above, one can obtain the calibration curves (E MPB vs |C| 2 ) for different cavity mode-exciton detunings by using this two-oscillator model and varying cavity-exciton Rabi-splitting. The data for the detunings used in the experiment are summarised in Supplementary Figure 6. Finally, the total exciton density for each exciton-photon detuning is deduced using formula (2) from Supplementary Note 1.2 (there we assume that trion fraction T = 0, since at high excitation density the trion-photon strong coupling is quenched).

The strength of trion-polariton nonlinearity
For a given energy of MPB the value of the Rabi-splitting Ω T at each power density can be deduced using the calibration dependencies in Supplementary Figure 2. The effective strength of trion nonlinearity responsible for the Supplementary Figure 6: Calibration curves. ELG 00 vs |C| 2 for the MPB with a two coupled oscillator model, obtained by varying the cavity-exciton Rabi-splitting from 0 to 19 meV for the four experimental cavity-exciton detunings. Bold solid sections of the curves correspond to the experimentally observed ranges of MPB peak positions for the corresponding detunings.
quenching of strong coupling is defined as a rate of decrease of the Rabi-splitting Ω T with n tot 5 : Alternatively, in the first order approximation β eff T can related directly to the redshift of the MPB branch in Fig. 3b of the main text as follows where η T = δ 2 P−T + ( Ω T ) 2 /( Ω T /2) and δ P−T = +5.3 meV is the energy detuning between the trion energy level and the lower polariton branch arising from coupling between photon and neutral exciton only (see above) Supplementary Figure 7: Nonlinear coefficient for trions. β eff T as a function of ntot. The error bars are 95% CI deduced taking into account the random error in the determination of the trion-polariton peak position at each power from the fitting procedure.
for n tot < 2 · 10 2 µm −2 . In Supplementary Equation (5) the minus sign explicitly accounts that Rabi frequency is decreasing function of density.
In the first order approximation our theory predicts a constant value of β eff T with density, or, in other words, a linear reduction of the trion-polariton Rabi-splitting with density (see Supplementary Note 3). As it is seen in Fig. 3c of the main text there is a good qualitative agreement between the dependence of the trion Rabi splitting and the theoretical prediction as a function of trion density: the experimental average value of β eff T = 37 ± 3 µeVµm 2 is in quantitative agreement with the theoretical estimate of 30 µeVµm 2 . Nevertheless, it is seen that the experimental points in Fig. 3c of the main text are not precisely positioned on the straight theoretical line. This is reflected by the fact that the experimental values of β eff T (deduced from the experimental data in Fig. 3c of the main text) vary from 120 to 20 µeV·µm 2 over the density range n tot from 0 to 200 µm −2 , which is shown in Supplementary Figure 7. We believe this variation of β eff T with n tot observed in the experiment may arise from the higher order effects due to composite nature of trions, not included in the theoretical treatment (see Supplementary Note 3).

Parameters characterising neutral exciton-polariton nonlinearity
The energy blueshift of the neutral exciton polariton arises from two mechanism: (1) the reduction of exciton-photon Rabi-splitting Ω X , and/or (2) blueshift of the neutral exciton level E X due to Coulomb exchange interactions.
Mechanism (1) is characterized by the rate of reduction of Rabi-splitting Ω X with exciton density 5 : It can be deduced from the MPB polariton energy blueshift assuming that only mechanism (1) is dominant: where η X = δ 2 C−X + ( Ω X ) 2 /( Ω X /2) and δ C−X is the detuning between bare photon mode and exciton level for a given data set.
Mechanism (2) is characterised by the rate of the blueshift of the exciton level with exciton density 5 : Similarly, assuming that only mechanism (2) contributes to polariton blueshift g eff X can be related to the polariton energy shift 5 : where ξ X = [1/2 + δ C−X /(2 δ 2 C−X + ( Ω X ) 2 )] −1 is the inverse of an excitonic fraction. In the experiment we can measure only the nonlinear behaviour (blueshift) of MPB. The UPB cannot be measured due to tunability of our laser. Therefore, experimentally we cannot separate the contributions to the neutral excitonpolariton optical nonlinearity from mechanisms (1) and (2). However, assuming that either only mechanism (1) or (2) is responsible for the blueshift of MPB, we can deduce the dependencies of the upper limits of β eff X = −d( Ω X )/dn tot and g eff X = dE X /dn tot factors on exciton density 5 , respectively. In the main text ( Fig. 5) we show that the theoretical g th X parameter is in semi-quantitative agreement with the experimental values of the upper limit of g eff X in the range of exciton densities 3 · 10 3 < n tot < 3 · 10 4 µm −2 . Now let us assume instead that the mechanism (1) is the only dominant mechanism over the whole density range. In this case we can observe that the theoretical β th X -factor is well below the experimental values of the upper limit of β eff X at n tot < 10 4 µm −2 as shown in Supplementary Figure 8. β th X approaches the experimental values only at higher densities n tot > 3 · 10 4 µm −2 . Such a discrepancy between the experiment and theory indicates that our assumption is incorrect; phase space filling for neutral exciton-polaritons [mechanism (1)] becomes important only at very high exciton densities n tot > 3 · 10 4 µm −2 , when the average distance between excited excitons is less than 5-6 nm and becomes comparable to the exciton Bohr radius a B ∼ 1 nm. By contrast, mechanism (2) is the dominant mechanism at intermediate exciton densities n tot < 3 · 10 4 µm −2 .
6. Measurement of the cavity mode decay rate (photon lifetime) The bare cavity mode decay rate (γ c ) was deduced from the temporal decay of emission intensity of the bare LG 00 mode (without flake in the cavity) excited resonantly with 100-fs pulse. The measurements were performed using streak-camera with the resolution time of 2 ps (see Supplementary Figure 9). We obtained the experimental cavity lifetime of approximately 4 ps, which corresponds to the FWHM of bare cavity mode ∼ 400 µeV, the value we measured on the spectrometer. a. Nonlinear refractive index n 2 of hybrid microcavity-MoSe 2 polariton system.
We note that applications of 2D materials imply that they would be integrated into photonic structures made of bulk semiconductors/dielectrics. Therefore, when characterising nonlinear optical properties of 2D materials it is useful to consider the nonlinearity of the whole hybrid 2D materials-semiconductor/dielectric photonic system and compare it with that of bulk photonic materials.
In order to compare trion-polariton nonlinearity with Kerr-like optical nonlinearity observed in bulk materials we can treat our open-access microcavity system with embedded MoSe 2 as a microcavity filled with a bulk of some nonlinear optical material, characterised by the effective nonlinear refractive coefficient n 2(MC) . The n 2(MC) coefficient due to trion-polariton nonlinearity can be estimated taking into account the effective optical path covered by a photon during the round trip between the two mirrors of the microcavity, which must be equal to an integer number of polariton wavelengths, Here n eff is the effective refractive index of the microcavity, N pol is the total number of polaritons excited inside the MC with a single pulse, L is the effective cavity length, λ is the wavelength of the trion-polariton emission in free space, δλ is the nonlinear shift of trion-polariton resonance, m is the order of the longitudinal cavity mode coupled with the trion. 2n 2(MC) Lλ L is the nonlinear optical path acquired by photon during the round trip between the two mirrors. Given absorption is the dominant process in our system N pol ≈ n tot .
Using Eqs. (5) and (10) we get the following expression where η T = δ 2 P−T + ( Ω T ) 2 /( Ω T /2), and δ P−T is the energy detuning between the trion energy level and the lower polariton branch arising from coupling between photon and neutral exciton only (see above). Taking into account the effective cavity refractive index n eff ≈ 1 (since most of the cavity electromagnetic field is confined in the gap between the two mirror), the wavelength of trion-polariton resonance λ ≈ 760 nm, the effective cavity size L ≈ 1 µm, β eff T = 37 µeVµm 2 and η T = 2 (δ P−T = 0) we estimate n 2(MC) ∼ 1.4 · 10 −13 m 2 W −1 . This n 2(MC) is about four to five orders of magnitude larger than 1.82 · 10 −17 m 2 W −1 in planar AlGaAs waveguides in the weak coupling regime 6 and 6 · 10 18 m 2 W −1 in silicon 7 and InGaP 8 , which have been used in a suspended membrane photonic crystal geometry. Kerr nonlinear effects (optical bistability) have been investigated in slab photonic crystal Si microcavities with embedded graphene layer 9 . The effective n 2 of hybrid graphene-Si microcavity system has been derived to be of the order n 2 ≈ 7.7 · 10 −17 m 2 W −1 , which is ∼ 3 − 4 orders of magnitude less than n 2(MC) due to trion-polariton nonlinearity. Finally, we note that the value n 2(MC) ∼ 1.4 · 10 −13 m 2 W −1 due to trion polariton nonlinearity is an order of magnitude higher than n 2 ∼ 1 · 10 −14 m 2 W −1 reported by us in neutral exciton-polariton GaAs-based system 10 .
To the best of our knowledge no Kerr-like nonlinear optical effects were studied in microcavities with embedded TMDC materials in the weak light-matter coupling regime. However, there were several studies of the effects associated with Kerr-like optical nonlinearity of bare layers of TMDCs (MoSe 2 , MoS 2 , MoTe 2 ) and graphene in the weak light-matter coupling regime on a picosecond timescale 11 . The values of n 2 coefficients for TMDCs layers were measured in the range 10 −16 -10 −17 m 2 W −1 depending on the excitation energy (above or below band gap).
The reference [12] reports the n 2 coefficient for WS 2 monolayer to be about 1.1 · 10 −15 m 2 W −1 . The value of n 2 coefficient for pure graphene flakes was measured about 2 · 10 −15 m 2 W −1 . 11 In Ref. [9] authors studied nonlinear hybrid Si-graphene microcavity, and deduced n 2 coefficient for a single graphene to be of the order 10 −13 m 2 W −1 .
In our trion-polariton microcavity system we can derive the effective nonlinear refractive index n 2(MoSe2) per single TMDC monolayer taking into account that in reality the nonlinear optical phase is acquired by light only on the passage of the monolayer during the round trip between the mirrors. Therefore, n 2(MoSe2) can be simply obtained by normalising n 2(MC) to d MoSe2 /L, where d MoSe2 ∼ 1 nm is the MoSe 2 thickness, yielding n 2(MoSe2) ∼ 1.4·10 −10 m 2 W −1 . This value is at least five (three) orders of magnitude larger than in TMDC 2D materials (graphene) studied in the weak light-matter coupling regime without formation of polaritons.

Supplementary Note 2
Nonlinear neutral exciton-polaritons: theory In the first section we describe the exciton-photon coupled system, accounting for the composite electron-hole (e-h) nature of neutral exciton. The case of a trion mode coupled to the optical mode is considered in the next section.

Exciton-polariton Rabi-splitting
To start, we consider an optical cavity described by the bosonic annihilation and creation operatorsĉ andĉ † , such that their commutation relations are [ĉ, Coupled to a semiconducting medium, an optical photon creates an exciton, corresponding to the bound electron-hole pair. The creation of an electron with the momentum q is described by the fermionic operatorâ † q , and hole creation is described by the operatorb † q . Their corresponding anti-commutation relations read where δ q,q is Kronecker delta function (and the same holds forb q ). Accounting for the attractive Coulomb interaction between an electron and a hole, the excitonic operator can be written as a composite bosonX ν , where ν is a general index which denotes the center-of-mass (CM) and internal degrees of freedom. The transition between electron-hole and composite exciton picture follows asX Here, k α,β correspond to electron and hole momenta, and k β , k α |ν is an exciton wave function written in the momentum space. The reversed transformation for describing an electron-hole pair in bosonic language can be written asâ † kαb † k β = ν ν|k β , k α X † ν , where summation goes over possible states of composite excitons ν in the appropriate orthonormal basis, such that ν |ν ν| = 1.
In the following we are interested in the system with strong light-matter coupling, where a composite exciton of certain CM momentum is coupled to the cavity mode. The Hamiltonian for the considered system readŝ where the first and second terms describe the free energy for the cavity photon mode,Ĥ cav = q ω cav,qĉ † qĉq ( = 1 hereafter), and composite excitonĤ X (i.e. coupled electron-hole) Hamiltonians. ω cav,q denotes the two-dimensional dispersion for the planar cavity mode, with typically ultralow mass, such that only small q's are considered. The third term describes the coupling between light and matter excitations. It can be written as a creation of an electron-hole pair by the cavity field with a coupling constant g, where in the second equality we have exploited excitonic form for the electron-hole pair, and considered dipolar transition with negligible transferred cavity momentum,ĉ q→0 ≡ĉ, being a usual assumption for description of strong coupling. Here, g = epcv m 2 2 0ωcav LcavA , where p cv is a matrix element for valence-to-conduction band transition, m is a free electron mass, L cav is a cavity length, A is an area of the system. We consider an exciton mode at the fixed center-of-mass momentum, which for brevity is set to zero,X 0 , and derive the corresponding Heisenberg equations of motion. It reads The first term generally describes the energy ω X at which the excitonic mode oscillates. The second, being proportional toĉ operator, provides the coupling to photonic mode, which we generally denote as G = g kα,k β i i|k β , k α X 0 ,X † i . If exciton corresponds to an ideal boson, i.e. [X,X † ] = 1, the coupling term reduces to G = g kα,k β i|k β , k α ≡ g k φ * k =: Ω X (0)/2, where we introduced the relative electron-hole momentum k and the Fourier transform of the exciton wave function φ k . This energy corresponds directly to the Rabi energy for the light-matter coupled system. Performing the diagonalization of the system at zero detuning (ω cav = ω X ), the Rabi-splitting between normal modes of the system is equal to Ω X (0).
We proceed by considering the composite structure of exciton, which is formed by two fermions. This comes from the fact that creation of a (correlated) electron-hole pair is not equivalent to a boson creation, as long as number of created pairs grows. It originates from the Pauli exclusion principle, which does not allow certain pair configurations in the full fermionic treatment, while disregarded in the purely bosonic picture 13,14 . The details for the difference between two cases were worked out by Monique Combescot and co-workers, and summarized in the so-called coboson approach to excitonic systems 13 . In the following, we apply use the coboson formalism to find corrections to Rabi and exciton energy appearing due to effects of non-bosonicity.
The main consequence of the composite nature of exciton is its peculiar statistics, which resembles bosonic one for small e-h pair concentration n ≡ N/A, but changes once it becomes comparable to the inverse of the effective exciton area. The generic commutation relations between composite bosons can be formulated as (see Ref. [13], eq. [4.16]) where an operatorD mi describes the deviation from bosonicity for excitons due to its composite nature.
In particular, this can be observed when one writes the commutator in Supplementary Equation (15) using the expression for composite exciton with zero CM momentum and relative momentum k, which is described byX The explicit form for the deviation operator isD 00 = k |φ k | 2 (â † kâ k +b † kb k ). Its structure thus hints that the deviation depends on the electron (or exciton) number N .
To calculate the influence of the non-bosonicity on Rabi energy renormalization we need to estimate the expectation value of the last term in Supplementary Equation (14) considering the (unnormalized) many-coboson state |N = (X † 0 ) N |ø , and singling out the prefactor in front of the cavity photon operatorĉ. Here, |ø denotes coboson vacuum state, and corresponding norm reads N |N = ø|X N 0 (X † 0 ) N |ø 1/2 . Note that in the case of composite bosons it was shown to differ exponentially from an ideal boson normalization for large N , 13 though for physical observables the difference appears as higher order terms in small N expansion.
The expectation value for the commutator can be written as The first term in Supplementary Equation (17) yields g k φ * k and is simply a Rabi frequency in the dilute system limit. The second term, however, involves the non-bosonicity operator. Accounting for the many-coboson state explicitly, it can be rewritten as where we have accounted for the fact that the action of the deviation operator on the ground state gives 0, i.e. D 0i |ø = 0 · |ø . Thus, its estimation relies on the commutator of the deviation operator with an exciton creation operator to the power N . For the first power, this can be derived as 13 where λ denotes the Pauli scattering element for input indices (i, j) and output (n, m). It reads explicitly where φ i (r α1 , r β1 ) is a generic coboson wavefunction written in the real space representation. It can be rewritten as a product of CM and relative motion component, φ i (r α , r β ) = (e iQi·R αβ / √ A) r αβ |i , where R αβ and r αβ are CM and relative coordinates for an e-h pair, which correspond to coboson with CM momentum Q i and relative motion quantum number i. The relative motion wavefunction can be also rewritten in the momentum space as r|i = k r|k k|i = k (e ik·r / √ A) k|i , which we will use in future. Proceeding with the estimation for the influence of the deviation in the many-coboson state, the commutator with N -exciton creation operator reads Using this, Supplementary Equation (17) can be rewritten as One can immediately see in the second term on the RHS that the Rabi frequency depends on exciton concentration N , times the prefactors coming from Pauli scattering elements. First, let us consider the intuitively easy case where the internal coboson index i coincides with the mode of interest, labeled as 0. Later, we show that this corresponds to the lower-order-in-N correction. In this case the expectation value ø|X N Considering the first Pauli scattering term λ(0, i; 0, 0), the summation over internal coboson index is performed as Here, for passing through the first equation sign we used that: 1) coboson states form a full orthonormal basis, i |i i| = 1; 2) transition element between between real and momentum space reads r α2 −r β2 |k = e ik·(rα 2 −r β 2 ) / √ A. For the second equality we exploited Dirac delta function definition in 2D, being dre ik·r = Aδ k,0 , which reduces summation to a single index, and recall our definition k|0 ≡ φ k . The second Pauli scattering term gives the same contribution.
In the case of 1s neutral exciton in a two-dimensional material the relation motion part of wavefunction in real space can be written as φ(r) = 2/πa 2 B exp(−r/a B ), with r being the relative (e-h) coordinate. The momentum space version then reads where a B corresponds to the 2D variational parameter.
Collecting everything together and performing summation as k → A (2π) 2 dk, the renormalized Rabi frequency as a function of concentration (in lowest order of na 2 B ) reads where G(0) = g 2A/πa 2 B = epcv maB 2 π 0ωcav Lcav . Finally, we need to account for the fact that G(n) is a derived term in the equations of motion, which includes ∝ n correction. Thus it originates from the effective nonlinear Hamiltonian H (nonlin) coupl which contains three excitonic operators and a photonic one, and on the contrary to the linear case provides extra factor of 2 in the equations of motion. The modified Hamiltonian with exciton density-dependent Rabi frequency then reads (at fixed momentum)Ĥ where renormalized Rabi frequency is and we denoted the exciton concentration as n X . The result above coincides with the estimates by Tassone and Yamamoto, 15 and Rochat et al., 16 although derived in a different way, without involving Usui transformation. We can proceed to calculate the terms being higher order in (na 2 B ). This relies on the exact calculation of average as 13 where F N is a coefficient which defines the deviation of statistics through N |N ≡ N !F N , with F N being 1 for purely bosonic states. Exploiting coboson theory, the ratio reads and The first term in Supplementary Equation (29) corresponds to the previously obtained case with ∼ na 2 B scaling. Performing the same procedure as before, we extend the results to include n 2 a 4 B contribution. After some algebra, we obtain where the derivation was performed up to ∼ n 3 a 6 B terms. The corresponding modified Rabi frequency then reads Interestingly, we observe that the light-matter coupling strength is a monotonically decreasing function of density, and its expansion has a sign-changing pattern. Thus, going in higher orders of nonlinearity shall help to obtain monotonically decreasing function.

Nonlinear exciton energy shift
In the previous section, we have considered the renormalization of the light-matter coupling term, which gains n-dependence. Now, let us consider similar effects which provide nonlinear energy term for composite excitons. To do so, we write the exciton Hamiltonian in terms of basic constituents, being electrons and holes. This is given bŷ where the first two terms correspond to energies of the electrons and holes, and readĤ = E g corresponds to the bandgap energy, and m e,h are effective electron and hole masses, respectively, measured in units of the free electron mass. The third term in Supplementary Equation (32) corresponds to the electron-hole Coulomb interaction, which ultimately leads to the formation of the bound state. It readsĤ e−h = − p,p ,q V qâ † p+qb † p −qb p â p , where V q = 2πe 2 /(AqS q ) is the standard Fourier transform for Coulomb interaction in 2D, S q denotes the screening function (so far unspecified), and we accounted explicitly for attraction between an electron and a hole. Finally, the termsĤ e−e andĤ h−h correspond to Coulomb interaction with only electrons and only holes. When accounting for the excitonic structure of e-h complexes, these lead to Kerr-type exciton-exciton interaction.
To account for the non-bosonicity and exchange-based Coulomb scattering between the cobosons on the equal footing, we can derive the total energy of the system, E X = Ĥ X , where expectation value is taken over N -exciton state. Following Ref. [13], the nonlinear (quadratic and higher) contribution to the energy of N excitons reads is a density independent exciton energy, and nonbosonicity factor reads is a direct Coulomb scattering term between excitons [see Supplementary Equation (20) for comparison], and V f f (r f1 , r f 2 ) corresponds to the real space Coulomb potential between carriers f, f . The ξ in term denotes exchange Coulomb scatterings, where either electron or hole is swapped between two composite excitons, which combines Pauli scattering and Coulomb interaction. Finally, the exchange term between three composite excitons (where carriers are swapped but no Coulomb vertex is included) yields For composite excitons with V αβ (r) = −V αα (r), which is true for the electron-hole potential, the direct term vanishes, ξ 0 0 0 0 = 0. This can be seen as well-known absence of direct contribution at zero exchanged momentum, valid both for III-V semiconductors and TMDs. Similarly, n λ 0 n 0 0 ξ n 0 0 0 = 0, and only exchange terms shall be accounted. They can be calculated as and The sums (37) and (38) can be converted into intergrals, and evaluated numerically for the exciton wavefunction in the form (24). As an important consequence of the monolayer structure of the TMD, the potential is chosen to be screened, [17][18][19] and has the form where e is an electron charge, 0 is a vacuum permittivity (note that SI units are used), r 0 is a screening length, and κ = ( s1 + s2 )/2 is an average dielectric permittivity for substrates from two sides. 19 This directly follows from the well-known Keldysh potential of the form defined with the help of Struve and Bessel functions of the second kind, and shown for the case of two electrons. Finally, collecting the terms up to N 2 order, the nonlinear energy of the excitonic mode (as appearing in the Hamiltonian) can be written as where exchange integrals for dimensionless length (x = r/a B ) depend on the screening length r 0 and read where in the case of TMDC materials the dimensionless parameter r 0 /(κa B ) enters the integrals. We perform the calculations considering MoSe 2 on hBN, where r 0 = 4 nm 18 , κ = ( s1 + s2 )/2 = 4 for hBN substrates, and leaving a B as a tuning parameter.

Exciton-polaritons at increasing density
Taking our previously derived results for the renormalization of coupling and exciton properties, let us translate it to the case of polaritons. 5 In the cases where trion mode can be excluded out considerations (it is weakly coupled and/or largely detuning), we use the two coupled modes Hamiltonian, concentrating on the exciton-photon coupling. This is justified by the experimental data at large detuning, where the trion fraction at X-C anti-crossing is estimated to be small (< 2%). The normal modes of neutral exciton-polariton system read where E − (n X ) ≡ E MPB (n X ) corresponds to the middle polariton mode we are interested in. We proceed with applying the presented theory to explain the nonlinear blueshift of the middle polariton branch as function of an exciton concentration, for the case where trion resonance is largely detuned. Using the coefficients derived above, and corresponding density dependences for the coupling Ω X (n X ) and exciton energy E X (n) terms, we plot the nonlinear energy shift as function of concentration. The theoretical results are shown in Fig. 4(b,d) of main text by red solid curves. Taking the exciton Bohr radius as the only fitting parameter, we set it to be a B = 0.85 nm, which allows to qualitatively the behavior of the system. Furthermore, we verify the obtained value performing the variational procedure to obtain exciton properties in MoSe 2 covered with hBN. This can be done using the standard procedure with screened Keldysh potential with r 0 = 4 nm, κ = 4, and the effective electron mass m e = 0.8m, 20 which was measured to be rather large in MoSe 2 and similar to the effective mass of the hole, taken m h = 0.84m. This leads to the reduced mass µ = 0.41m. Performing minimization, we get a = 259 meV. These values lie close to the fitted value and experimentally measured energy, respectively. (As a bonus, in the next section devoted to trions we explicitly show how the same result can be easily obtained in the momentum space. ) We note that the ability to reproduce measured energy shift of E MPB is only possible once both Rabi frequency reduction and nonlinear interactions are considered, while otherwise failing to provide required scaling. Namely, the inclusion of Coulomb-based exchange can only explain the observed behavior (2 meV shift within at order of magnitude change for the density) for either largely increased exciton Bohr radius in TMDC, which is unlikely, or much higher concentration going into ∼ 10 13 cm −2 range. At the same time, if only Rabi renormalization is accounted, the saturation of nonlinear shift cannot be reproduced.

Trion-polariton Rabi-splitting
We consider the system with an initial doping, and study the effects of light-matter coupling with a multiparticle bound state. In the MoSe 2 TMD this corresponds to a negatively charged exciton -trion -which is spectroscopically located 30 meV below the excitonic resonance. We aim to estimate of the Rabi-splitting change for the case of a trion. In the similar fashion, the deviation from ideal statistics changes the value of the trion-photon coupling. However, we note that the strong light-matter coupling regime for trion is much less studied, and its treatment requires extra care.
We begin with the interaction between the cavity and the trion mode. The latter can be generally described by a composite creation operator which creates two electrons and a hole from the vacuum state,â † ke,seâ † k e ,s e b † k h ,s h |ø . Here, k e,e ,h are the momenta of the respective individual constituents (so-called carrier coordinates 21 ), and s e,e ,h are the spin indices. We note that the most favorable trion configuration in MoSe 2 monolayer is the singlet state with two electrons having anti-parallel spin. 22 The operator corresponding to the creation of a singlet state can be written as 23T † with wavefunction being separated into a trivial center-of-mass part with momentum K, β e = 1−β X = m e /(2m e +m h ), and the relative motion part described by trial wavefunction for the relative motion φ T k1,k2 (we consider zero CM momentum case). The wavefunction is written for the relative electron-hole coordinates r 1 and r 2 , being radiusvectors between the first electron and the hole, and the second electron and the hole, respectively. In the real space it corresponds to the two exponentially decaying functions where λ 1 and λ 2 are the variational parameters corresponding to the distances between electrons and a hole. Note, that φ T (r 1 , r 2 ) is symmetrized and is normalized to unity, with χ = 4λ 1 λ 2 /(λ 1 + λ 2 ) 2 . The momentum space version then reads where we defined N = [2(1 + χ 2 )] −1/2 and φ (j) Finally, reordering the electron operators in Supplementary Equation (45), and considering CM momentum much smaller than typical relative momenta, we can write trion creation operator aŝ The choice of the wavefunction (47) is of course far from optimal, as to describe quantitatively the shape of trion solution more complicated ansatzes with hundreds of orbitals shall be used. 24 However, in order to get any sensible result for Rabi frequency renormalization, this is the form we shall adopt.
Once there is a non-zero number of free electrons, the absorption of a circularly polarized photon can then allow a creation of a trion. The Hamiltonian of the system can be written as the sumĤ T =Ĥ T 0 +Ĥ T coupl of non-interaction cavity/electron/trion HamiltonianĤ T 0 and the coupling Hamiltonian for light and matter, with g being conduction-to-valence band transition matrix element, previously defined in the exciton case. ε k is an electron dispersion, ω T K is a trion dispersion, and in Supplementary Equation (49) the summation over spin is assumed. In Supplementary Equation (50), similarly to the exciton case, the wavefunction of the relative motion appears due to the fact that out of free electron-hole complex the bound trion state appears. The process in Supplementary Equation (50) in simple terms can be seen as a creation of the electron-hole pair attached to an electron in the Fermi sea, while the electron state is (slightly) changed. This can be conveniently described by the quasi-bosonic excitation, defined by the operatorB j , which readsB † It creates an excitation out of Fermi sea state |FS , and N e is a number of free electrons available for the trion creation, which can correspond to the selected spin configuration, meaning that the total number of electrons in the system reads as N tot e = 2N e . The combinatorial prefactor 1/ √ N e comes from the number of different ways the excitation can be created, 25 and we note that (51) holds for low temperatures where electron gas is degenerate. As the excitation operatorB † K represents a composite boson, similarly to excitons described in the previous Supplementary Note 2, it is prone to the phase space filling effects. At the same time, it is not a bound state, and thus exhibits different statistics deviation behaviour. We will describe this point in details later, when trion-based saturation effects are considered.
We note that previously the light-matter coupling in MoSe 2 TMD material was also considered for the case of exciton-polarons. 26 This corresponds to similar creation operatorB † j , but different ansatz for the wave function, which accounts for dressing of photo-created exciton with electrons in the Fermi sea. We stress that both triondominated or polaron-dominated regimes can be possible, as considered in Ref. [27]. Recent predictions for the GaAs samples estimate the cross-over into exciton-polaron regime to happen for Fermi wave vector to be comparable with inverse Bohr radius for the system, k cr ∼ 0.8a −1 B . 28 As we show later, the experiment is conducted in the k F a −1 B limit of small concentration, which corresponds to the trion-dominated regime. Next, we proceed with the estimation of trion Rabi-splitting, or conversely the free electron density. It is given by the bare coupling constant g multiplied by the square root of electron density and the wave function part responsible for absorption renormalization due to the confinement. First, we account for the fact that photon momentum has typically small values q, being much less than other relevant wavevectors. This also translates into nearly zero centerof-mass momentum of the generated exciton, where an electron and a hole are located close to each other. Here, we follow the approach introduced in Refs. [29,30], where the so-called electron-exciton coordinates are used. These correspond to the electron-hole relative coordinate (seen as an exciton) and the relative coordinate of the CM for "exciton" to the second electron. They are generally described by length parameters λ and λ . Due to complex symmetrization requirements, they do not allow to choose wavefunction in the simple form, thus preventing the analytical calculation. However, in the limit of large λ 2 λ 1 the e-h and e-X coordinates become nearly equivalent, and we can set λ ≈ λ 1 and λ ≈ λ 2 .
Taking the trion wavefunction to be Fourier transformed with respect to an exciton internal motion, the coupling then can be estimated setting r 1 = 0 (i.e. for the closely located photocreated e-h pair). Simultaneously, we shall account that Fermi wavevector k F λ −1 1 , λ −1 2 , and thus the wavefunction in Supplementary Equation (50) can be considered as a constant at k 2 λ 2,1 = 0, going in front of the sum. Altogether, the coupling Hamiltonian can be rewritten asĤ where the trion Rabi frequency reads where we used the Rabi frequency definition for the neutral exciton case, Ω X /2 = g 2A/πa 2 B , and n e = N e /A is a concentration of free electrons which can form trions. Supplementary Equation (54) then allows to estimate n e using once the variational parameters are known.

Trion binding energy and variation
Next, we proceed to define λ 1 and λ 2 for trions in TMDC. We follow the approach outlined in Ref. [23], using the wavefunction (48). The expectation value for the trion Hamiltonian (includes kinetic terms for relative motion and Coulomb interaction) then can be written as where we define auxiliary quantities: where γ = m e /m h ,λ = λ 1 λ 2 /(λ 1 + λ 2 ), and we remind that r 0 is screening parameter, κ is average dielectric permittivity of the substrate. Here, all length parameters are measured in the units of and energies are measured in units of To obtain the binding energy for the trion complex, we minimize E (T) [λ 1 , λ 2 ] with respect to variational parameters, and subtract the exciton binding energy contribution E X b . The latter is obtained from minimization of which for previously defined parameters of m e = 0.8m, m h = 0.84m, r 0 = 4 nm, κ = 4, gives λ 0 = 0.93 nm and binding energy of E X b = − min{E (X) [λ 0 ]} = 259 meV. The variational procedure for the trion then gives the binding energy of E T b = − min{E (T) [λ 1 , λ 2 ]} − E X b = 26 meV for λ 1 = 0.87 nm and λ 2 = 2.54 nm. These are the parameters which will be used in the following. Although we remind that considered variation with two parameters is oversimplistic, it provides energy estimate to be very close experimentally measured trion binding energy of 30 meV.
Finally, substituting obtained radii λ 1,2 , a B = 0.93 nm, and experimentally measured Rabi-splittings Ω T = 5.8 meV, Ω X = 17.2 meV, using Supplementary Equation (55) we estimate the electron concentration available for trion creation to be n e = 4.05 × 10 10 cm −2 , with the full concentration corresponding to n tot e = 8.1 × 10 10 cm −2 . Note that so far only variation and parameters obtained from the exciton-polariton case were used, with no fitting involved.

Trion Rabi-splitting quench
We continue with the calculation of the modified trion Rabi frequency due to the deviation of statistics. For this, similarly to excitonic case [Supplementary Equation (14)], we derive the equations of motion for the excitation modê B j using Hamiltonian (53). The nontrivial dynamics part comes from the light-matter coupling term and relies on the calculation of commutator [B q ,B † q ]. To do so, it is instructive to rewrite the excitation operator in terms of trion and electron, yielding In the case of low pumping intensity, the trion anti-commutation relations resemble that of ideal fermions, {T i ,T † j } = δ i,j , and the number of trions goes to zero. Then, the commutator for q = q reduces to integral over distribution function f k and gives unity, as it should be for an ideal bosonic mode. However, for the increase of pumping we observe two contributions which change the commutation relation. First contribution comes from deviation of fermionicity for composite trion operator, such that {T q +k ,T † q+k } is not a simple delta function anymore. For equal momenta this starts at the value of unity, and decreases with powers of n T λ 2 .
The second contribution comes from the composite nature of the quasi-bosonic operatorB q , which ultimately depends on the ratio between number of trions and available free electrons for their creation. To demonstrate this point, let us rewrite Supplementary Equation (65) as where the deviation operator is formally introduced aŝ where we have considered the limit of n T λ 2 1, such that trions can be approximately treated as ideal fermions. We observe that while conceptually the deviation operator resembles the one used for excitons in the previous section, the fact that trion-electron excitation is not bound leads to different closure relations and statistics. In this case, it is reminiscent to intersubband excitations 32 with trion being an excitation over the Fermi sea. We proceed by deriving the commutation relations for the devation operator and excitation operator, which reads and it can be recursively generalized to the case of N T particles as The derived commutation relations, which are dependent on N T /N e ratio, will be later shown to ultimately lead to the quench of the trion Rabi frequency, where the commutator B q ,B † q averaged over highly-excited many-body state vanishes.
In the following, we proceed with considering the two aforementioned contributions one-by-one.

Deviation from fermionicity for the composite trion anti-commutator
To account for the deviation of statistics, we use the generalized many-body formalism for composite n-particles 31 . This allows to calculate anti-commutator for the composite fermion (trion in our case), which consists of three particles of different flavor (opposite spin electrons an a hole). In the general form, it reads where deviation from fermionicity operatorΞ mi is defined as and the operatorΞ † mij is defined through the anticommutator Here single exchange integrals λ ρ (j, i, n, m) are as in Supplementary Equation (20) (though with three particle wavefunction), and ρ denotes the carrier to be exchanged. In total, it provides six contributions, which in the case of zero exchanged momentum we expect to be the same. The last term in (72) corresponds to three-particle exchanges Λ p,k,n,j,m,i with all permutations. As it typically corresponds to ∼ λ 4 scaling, which shall be accompanied by the quadratic density contribution, we refrain from considering it, and concentrate on lower order terms only. However, we note that the missing term might still effect the quench of the trion Rabi-splitting, as at increasing concentration the terms in all orders become important. The trion exchange is calculated using the trial wavefunction in the electron-exciton basis without symmetrization, which is accounted at latter stage, and we take hole exchange as an example. In analogy to the case of trion interaction 33 it reads To evaluate (75) we use new coordinates p − p = δp, (p + p )/2 = P, make momenta dimensionless multiplying it by λ, and define ξ = λ /λ. Finally, following the same procedure as for composite excitons, the anti-commutator N T |{T q ,T † q }|N T is averaged over a state of N T composite particles (now fermions), such that N T contributions is obtained. This leads to the estimate for the deviation where the dimensionless exchange integral reads and we account for symmetrization between λ and λ parameter as a separate term. As in the case of excitons, we observe that the first order correction in n Tλ 2 reduces the coupling as a function of the composite particle density (λ is an effective parameter of length dimensionality).

Effects of medium saturation
Next, we find that the increased number of trions as a consequence of increasing pump intensity provides another contribution for the trion Rabi frequency reduction. To derive the trion-density dependence for Ω T , we employ the same strategy as previously used for composite excitons in Supplementary Note 2.1. For this, we consider the ground state of the system as a Fermi sea of free electrons |FS available for the trion creation. The relevant excited states then correspond to multi-trion states |N T ≡ (B † q ) NT |FS . The nonlinear contribution to the trion Rabi frequency associated to the composite nature ofB q comes from the average where we used Supplementary Equation (70) and the fact that q, q , q are small. The analysis of Supplementary Equation (79) shows that modified commutation relations can ultimately leads to the quench of the strong coupling once the number of trions becomes comparable to half the number of free electrons N T = N e /2. However, this only corresponds to the lowest order corrections, and higher terms shall be accounted for increasing N T /N e ratio to get the full treatment. In this case, smooth reduction of Ω T is expected up to N T = N e . Finally, collecting all contributions together (including the one described in Supplementary Note 4), we can write the effective commutator at growing trion density n T = N T /A as a function where similarly to exciton case we conjectured the appearance of the quadratic term ∆ 2 nF /2, which appears in the expansion of the exponent. The function f T (N T , N e , λ, λ ) is decreasing from 1 to 0, and we consider it zero after the quench. The important parameter then is the half of available electron density n e /2, which defines the excitation density at which quench is observed. The resulting density dependent trion Rabi frequency, defined as in Supplementary Equation (64), then reads and is used later to calculate the density dependence for the polariton modes. Intuitive explanation. To explain the leading trend of linearly decreased coupling for the number of trion equal to half the number of free electrons, we note that the coupling of the trion mode to the cavity bares analogy to the atom-photon coupling 25,34 (contrary to the neutral exciton case). This is readily seen in the ∝ √ n e dependence for the coupling constant, similarly to the common square root enhancement for the N two-level emitters. 35 Given this correspondence, we show how the coupling between trion and photon changes for high excitation power. For the trion case, the state |g corresponds to a free electron, while excited state |e corresponds to the created trion. The total number of available excitations is thus equal to the number of free electrons n e ≡ N . Using the analogy, we write the Hamiltonian of the system asĤ whereâ † (â) is a creation (annihilation) operator for a cavity mode, and j corresponds to the considered two level system. The coupling term ∝ g thus describes polaritonic physics. For simplicity, we can take ∆ = 0 and measure cavity mode energy ω c from this value. The usual way to treat light-matter coupling in (82) is to assume weak excitation conditions and perform effective bosonization 35 (i.e. make Holstein-Primakoff transformation). For this, the excitation creation operator readŝ where 1/ √ N corresponds to the normalization condition, andb can be written similarly. The overall meaning ofb † is the creation of excitation out available two-level emitters (free electrons) as a superposition. With the new operators Hamiltonian (82) can be recast in the familiar form where the last term corresponds to the usual polaritonic coupling with the superradiant enhancement, as for coupling to an ensemble of emitters. To see the influence of the light-matter coupling on the energy of the system (polaritonic shift), we take the many-body wave function in the form where n ph corresponds to the number of photons and n exc to the number of excitations (i.e. number of |1 j atomic states). Only certain states will be coupled by the off-diagonal light-matter interaction term. The expectation value for the Hamiltonian (84) yields Ĥ = ω c n ph + n ph − 1, n exc |g √ Nb †â |n ph , n exc − 1 = ω c n ph + g √ N n ph − 1, n exc | √ n ph √ n exc |n ph − 1, n exc (86) where we considered number of excitations to be equal to number of photons. One observes that the expectation value contains the same value of coupling g √ N as before, which is equal to Rabi frequency. The situation is however different when weak excitation conditions are not met. In this case effective bosonic picture (and Holstein-Primakoff transform) is not longer valid. To compare the two cases, the full basis in (82) must be considered. Choosing the wave function with N − m states to excited, and m states be in ground state {0 m i }, we can write it as We note that (87) is of course an approximation, and the full state may contain components with different number of excited stated. However, as the coherent state distribution is expected, their influence is suppressed. Taking the expectation value for (82) for (87) one gets where the coupling reduced by the factor m/N . For growing number of excitations n exc (and increasing number of atom in ground state), this prefactor can be expanded into series, leading to ∝ m/N ≈ 1 − 2n exc /N + O[n exc /N ] 2 dependence. Finally, it is easy to see that for all states being excited (m = 0) the coupling goes to zero, as there are no states to couple. This leads to the conclusion that for the inverted medium (maximal number of trions), the cavity becomes decoupled from the matter. Of course, in reality other effects coming from non-fermionicity play the role, and linear quech to zero is expected to change into a smooth function.

Trion polaritons at increasing density
Finally, we describe the process of trion Rabi quench, evidenced by the behavior of the middle polariton branch. To describe Fig. 3(b) in the main text we perform the diagonalization of the full photon-trion-exciton system, which can be written as 2 Ω X (n X ) 1 2 Ω T (n T ) 1 2 Ω X (n X ) E X (n X ) 0 Here we consider the nonlinear contributions to both trion and exciton modes, taking experimentally measured values E (0) T = 1621.02 meV, E C = 1630.32 meV, E (0) X = 1645.72 meV, Ω T = 5.8 meV, Ω X = 17.2 meV, obtained exciton Bohr radius a B = 0.85 nm, variation parameters λ 1 ≈ λ = 0.87 nm, λ 2 ≈ λ = 2.54 nm, and using the estimated electron concentration n e = 4.05 × 10 10 cm −2 . The result is shown in Fig. 3b [main text] by the red solid curve, and allows to reproduce the initial red shift of E MPB due to quenched trion Rabi-splitting, followed by the weak blue shift caused by exciton-exciton interactions.
Finally, we notice that in experiments the energy blue shift of the MPB at small exciton concentrations (see Figs. 3b and 5 of the main text) is about one order of magnitude larger than that predicted by the theory. The missing blueshift contribution, which manifests as a plateau due to competition with Ω T quench and subsequent linear growth, is identified as a trion-exciton interaction. Indeed, the full trion Rabi quench appears when high order terms are neglected, with function f T (n T ) reaching zero non-smoothly. However, their account shall provide nonzero Ω T even for n T ≈ n e /2, leading to residual coupling and small trion admixture. This for instance can lead to several percent fraction of the trion in MPB branch, and contribute as trion-exciton exchange. The interaction can substantially increase the energy of MPB state even for the trion being weakly coupled to light. As a result, at exciton densities 5-10 times above the maximum density of the excited trions, the MPB can exhibit substantial blueshift well above that predicted by the theory, which considers only neutral exciton-exciton interactions.
The full calculation of trion-exciton interaction energy is formidable and is typically limited due to its large dependence on the trion wavefunction ansatz 33 . Nevertheless, the overall strength of trion-exciton interaction is estimated to be one-to-two orders of magnitude stronger than that of exciton-exciton exchange. This is because of the following reasons: 1) Firstly, the number of electron and hole exchange processes is increased compared to the X-X exchange case; 2) the outer shell of the trion, described by λ 2 , defines the scattering cross-section, which is larger than for exciton; 3) the direct Coulomb term is expected to play a role, unlike for the neutral exciton-exciton interaction case.

Supplementary Note 4 Quantum trion-polaritons: theory
We now present a theoretical estimate of how a system analogous to the one considered in the current paper can be used to observe strong nonlinear response at the quantum level of a few photons. The experimental signature of such behavior is a pronounced antibunching of the photon emission, which can be monitored by measuring the second order coherence function. To calculate the latter, we consider a system described by the trion-photon Hamiltonian, where the coupling is between a single photonic mode of the cavity, described by the bosonic operatorsĉ,ĉ † , and a trion mode with zero in-plane momentum K = 0, characterized by the operatorsB 0 ,B † 0 [see definition in Supplementary Equation (51)]. The Hamiltonian can be written in the rotating frame aŝ where P denotes the pump strength for a cw coherent optical drive of frequency ω p . In Supplementary Equation (90) the first two terms correspond to free cavity photons and trions and the third term describes the photon-trion coupling. The Rabi splitting Ω T is given by Supplementary Equation (54). The effective polariton energies in the weak excitation limit read ω L,U = (ω c +ω T )/2∓ Ω 2 T + (ω c − ω T ) 2 /2. Finally, it is convenient to go to the rotating frame with respect to the last term in Supplementary Equation (90), such that the system is described by the detuning from the pump frequency.
The non-bosonic operatorB † 0 is characterized by its matrix elements in the Hilbert space spanned by the Fock states for trions, |N T . The matrix elements read where F NT ≡ Ø|B NT 0B †NT 0 |Ø /N T ! corresponds to the correction factor accounting for the composite nature of the trion excitation (in the case of bosonic excitations F NT = 1). Its explicit form can be derived using the recursive relation which is similar to those reported earlier for Frenkel excitons 36 and intersubband excitations 32 . As for the numerator, it can be evaluated using the commutation relation derived iteratively from Supplementary Equation (65). A careful treatment of recursion in arbitrary order gives a closed expression for the matrix elements and a similar expression can be derived for its complex conjugate. Importantly, Supplementary Equation (94) works for the relevant case of few trion excitations N T ≤ N e (contrary to the Holstein-Primakoff approach 38-40 which fails in the limit where N T = 1). This will be important for obtaining correct quantum statistical properties for the cavity emission. The dynamics is studied by numerically solving the master equation for the full density matrix of the system in the truncated trion-photon Hilbert space. It reads where the first term on the right hand side corresponds to the coherent part of the evolution, the second term describes photonic dissipation (characterized by the finite broadening of the cavity mode related to the finite lifetime of the cavity photons), γ c = τ −1 c , and the third term describes trion dissipation characterized by nonradiative broadening γ T = τ −1 T . To characterize the statistics of the cavity output we evaluated the second order coherence function at zero delay for the cavity photons, defined as whereρ s denotes the steady state density matrix for the continuously driven dissipative system. In the current experiment we measured a Rabi frequency Ω T = 5.8 meV and estimated the free electron density n e = 4 × 10 10 cm −2 , with the corresponding number of electrons being ∼ 1200 in the cavity of area A = 3 µm 2 . These experimental values already imply the expected trion-photon coupling strength for the case of a single electron in the cavity area, g c = Ω T / √ N e = 0.17 meV. In improved devices the photonic mode area can potentially be reduced to A = 1 µm 2 using a curved top mirror with smaller radius of curvature 43,44 so that the coupling strength is enhanced by a factor of √ 3, giving g c = 0.29 meV. Meanwhile, the photon decay rates for an open cavity routinely reach ∼ 0.1 meV values, and can be as low as γ c = 10 µeV (∼ 65 ps lifetime) in the state-of-the-art samples 42,43 .
Furthermore, electron concentrations down to n e = 10 10 cm −2 can be realised in gated samples. In this case the expected Rabi splitting will be ∼ 3 meV for an average cavity occupation of 100 electrons, making phase-space filling effects even more pronounced and going far beyond the electronic confinement regime. Finally, in TMDC samples of high purity, which are encapsulated between thick hBN-layers, the trion inhomogeneous broadening as well as non-radiative recombination may become negligible such that the trion nonradiative linewidth γ T will be determined by pure dephasing due to scattering with phonons 41 . This may result in nonradiative trion linewidths as small as γ T ∼ 10 µeV at a temperature of 1 Kelvin 41 .
While studying the second-order coherence in the system, we note two possible mechanisms which can reduce the multi-photon component and facilitate single photon emission. The first mechanism corresponds to the conventional blockade-type antibunching, where two-photon occupation is suppressed by strong trion-photon coupling at the single particle level, g c /γ c,T 1. In the following we show that this shall be possible in future high-quality samples. The second mechanism can be identified as an unconventional-type single photon blockade 46,47 due to phase space filling effects, which does not require strong coupling and works at optimal parameters of g c ∼ γ T and ω p ≈ ω c (see 45 for the full analysis). It relies on destructive interference between the direct coherent optical excitation two photons and the trion-mediated excitation path, 47 thus relaxing the requirement for strong energy shift. At the same time, it causes oscillations of the second-order coherence as a function of delay, and generally has smaller emission probability. Below, we consider the two regimes as long-term and near-term goals for nonlinear trion-polaritonics with TMDC materials.
The results of the second-order coherence calculations are shown in Supplementary Figure 10. First, in Supplementary Figure 10(a) we plot g (2) (0) for a range of pump frequencies close to the cavity transition. For this, we consider cavity linewidth γ c = 0.05 meV, ω c = ω T , and g c = 0.29 meV with 100 electrons. Studying the dependence for different values of the nonradiative trion decay rate γ T we observe the appearance of an antibunching window when g c ≈ γ T . At the same time, for narrow linewidths γ T g c the antibunching behaviour disappears from the ω p ≈ ω c region, signifying the resonant interference-based nature of the effect and the modest coupling requirement. At the same time, we note the limited efficiency of the single photon emission in this window, as cavity occupation is typically in ĉ †ĉ ∼ 10 −3 ..10 −4 . We envisage that in practise the optimal parameters would be tuned using sample positioning in an open cavity as a tool.
In Supplementary Figure 10(b) we consider a different range of pump detunings, where the coherent drive is nearly resonant with the lower trion-polariton, ω p ≈ ω L , which corresponds to non-zero ω c − ω p detuning. In the calculations we assume a high-quality cavity with γ c = 10 µeV and improved values of the trion linewidth limited by thermal effects. The plot for g (2) (0) shows an antibunching for pump frequencies slightly below the transition, with single photon purity gradually improving as the coupling ratio g c /γ c,T increases. The Fano-shape profile of the g (2) (0) dependence draws the connection to conventional Kerr-based polariton blockade, 48 which can be accessed in structures of larger lateral size and favours the strong binding energy limit for excitons.
Finally, Supplementary Figure 10(c) shows the minimal value of g (2) (0), minimized over a wide range of pump detunings spanning both regimes, as a function of the trion nonradiative linewidth γ T . For trion decay rate comparable to the light-matter coupling constant g c the unconventional antibunching can be observed (γ T > 0.1 meV). For longlived trions with small non-radiative decay (0.01 meV) we also observe pronounced antibunching due to conventional blockade, which is limited by the cavity quality factor. This shows that single photon emission with trion-polaritons in MoSe 2 is possible in high quality samples.