Path integral approach to bosonisation and nonlinearities in exciton-polariton systems
Anna M. Grudinina Nina S. Voronova




National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Kashirskoe shosse 31, 115409 Moscow, Russia
Russian Quantum Center, Skolkovo IC, Bolshoy boulevard 30 bld. 1, 121205 Moscow, Russia


nsvoronova@mephi.ru
National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Kashirskoe shosse 31, 115409 Moscow, Russia
Russian Quantum Center, Skolkovo IC, Bolshoy boulevard 30 bld. 1, 121205 Moscow, Russia


Large exciton-polariton optical nonlinearities present a key mechanism for photonics-based communication, ultimately in the quantum regime. Enhanced nonlinear response from various materials hosting excitons and allowing for their strong coupling with light is therefore the topic of intense studies, both in theoretical and experimental domains. Reports on the scattering rates arising due to various system's nonlinearities, such as the exciton-exciton Coulomb interaction and the Pauli blocking that leads to the saturation of the exciton oscillator strength, however, are contradictory. In this work, we develop a formalism allowing to track the exciton nonlinearities appearing in the regime of strong coupling with photons, that includes finite temperatures, mixing of the exciton excited states, and the dark exciton contributions to saturation self-consistently. The equilibrium path integration approach employed here to address the polariton composite nature, leads to a transparent hierarchy of various contributions to nonlinearity. At the same time, by taking the simplest limit of zero temperature and so-called “rigid” excitons, through our framework we retrieve the expressions derived in conventional approaches for exciton interaction constants. In particular, our theory allows to clearly show that such interaction constants cannot be used as fitting parameters tunable in a wide range of values, as they are strictly defined by the material properties, and that other explanations are due for large optical nonlinearities recently reported.


Path integral approach to bosonisation and nonlinearities in exciton-polariton systems
    Nina S. Voronova
    ^1Division of Space Research and Planetary Sciences, University of Bern, Switzerland

    ^2Max Planck Institute for Astronomy, Heidelberg, Germany

    

==============================================================================================================================================================





§ INTRODUCTION




Exciton-polaritons, hybrid semiconductor quasiparticles that result from strong coupling of electronic excitations with light, are attractive candidates to endow photons with strong nonlinearity, both in the macroscopically coherent regimes and at the single-particle level <cit.>. 
This promising perspective, together with quickly developing state of the art, led in the last years to intensification of experimental and theoretical studies aimed at the analysis of polariton interactions in various material systems, including both conventional semiconductor quantum wells (QWs) <cit.> and two-dimensional (2D) transition-metal dichalcogenides (TMDs) <cit.>, the latter allowing to operate at elevated temperatures.


When considering macroscopic phenomena on the polariton level, such as the polariton Bose-Einstein condensation <cit.> and superfluidity <cit.>, excitons are most often treated as proper bosons, due to the underlying exciton-photon strong coupling. In the microscopic descriptions on the exciton level, however, it has been argued that their composite electron-hole structure needs to be taken into account. In this context, for electron-hole-photon systems with strong coupling, such as polaritons in a microcavity, several formalisms considering the effective exciton Bose field have been developed, e.g. directly introducing the exciton operators <cit.>, performing the Usui transformation <cit.>, or basing on the semiconductor Bloch equation <cit.>. All these methods represent generally the same approach describing the bosonisation of interacting fermionic constituents of the exciton in presence of the photon field in the case of low densities (n_ exa^2≪1, where n_ ex and a are the 2D exciton density and Bohr radius, respectively). On the other hand, the opposite limit where the coupling to photons is considered dominant and Coulomb electron-hole interaction absent, was studied within the generalized Dicke model in a series of works by Littlewood and co-authors <cit.> and, later, in Ref. <cit.> for excitons bound by light. In both of these limits, however, the considerations were mostly limited to zero temperature and spinless electrons and holes (except for Ref. <cit.>).

It is perhaps instructive to recall that for bulk semiconductors,
the electron-hole systems without strong coupling to photons were originally studied by Keldysh and Kozlov <cit.> in 1960s using the Green's functions technique with additional corrections from multiple electron-hole scattering processes. Later, Kiselev and Babichenko <cit.> proposed
the description of such electron-hole fermionic system based on the path integral approach. This theory allows for perturbative derivation of the effective exciton action for a weakly-interacting exciton Bose gas starting from the electron-hole formalism. The exciton interaction constant calculated in <cit.> and <cit.> is equal to g_ ex^ 3D = (13π/3) ħ^2 a_ 3D/M, which corresponds to 1s-exciton scattering (here M is the exciton mass).

In an analogous approach for quasi-2D excitons in semiconductor QWs, the exciton interaction constant g_ ex^ 2D = (1-315π^2/4096)4πħ^2/M was derived in 1995 <cit.>. Further works treating QW exciton interaction already in the polariton regime <cit.> reproduced this result and, furthermore, underlined another possible source of the polariton optical nonlinearity: the saturation of exciton oscillator strength. 

With the emergence of TMD monolayers as a versatile platform hosting excitons and exciton-polaritons, the search for correct estimates of the interaction constants g_ ex and g_ sat, characterising the exciton pair interaction and saturation, respectively, was renewed (here and below we drop the superscript “2D” for clarity). On the one hand, the small exciton Bohr radii and large binding energies hold promise for robust highly-nonlinear TMD polaritons up to room temperatures. On the other hand, the lowest-lying s-states of such excitons exhibit significant deviations from conventional Coulomb-bound hydrogenic model <cit.>, since the interaction between charges in atomically-thin layers is described by the Rytova-Keldysh potential <cit.>.
Among theoretical works, Ref. <cit.>,  revisiting the logarithmic factor in the interaction strength—a well-known result for 2D boson scattering problem <cit.>—underlined that for polaritons the interaction constants should be larger than those for excitons, due to the difference between the exciton and lower-polariton (LP) energies. Ref. <cit.> in their Supplemental Information provides the detailed derivation of the Rabi splitting renormalisation (saturation) in the operator formalism, similar to <cit.>, but up to the second order in n_ exa^2. Thus, the analytical expressions for the constants g_ ex and g_ sat describing both polariton nonlinearities, valid for the case of robust 1s–excitons at T=0, exist.

Nevetheless, recent experiments reveal not only a substantial disagreement with those predictions, but also contradictory results in comparison to each other. In particular, while some giant exciton saturation constants were recently reported for 1s-excitons in WS_2 <cit.>, Ref. <cit.> reports on the numbers two orders of magnitude smaller for the same material. Still, the saturation-related nonlinearity g_ sat in both works exceeds g_ ex (in Ref. <cit.> the nonlinearity due to exciton pair interaction is considered overall absent). Similar studies performed in <cit.> for MoSe_2 revealed g_ ex and g_ sat to be of the same order, whereas in <cit.>, these constants are smaller by one and two orders of magnitude, respectively. Notably, a very recent experiment performed with GaAs QW polaritons, which was able to resolve the upper-polariton (UP) branch, also revealed g_ sat≫ g_ ex <cit.>. In all the above experimental works the two constants are used as free parameters to fit the observed dispersions shifts with density, in disconnection with the theoretically-obtained numbers. All these results call for interpretation and detailed analysis.

In this work, we develop the path integral approach <cit.> for an electron-hole-photon mixture, examining the influence of the photon field on the shape of the exciton field, and derive the effective exciton-photon action. While the considerations presented here in their simplest limit (T=0 and the rigid 1s-exciton) do not lead to results other than those obtained in previous works <cit.>, they differ in regard to their derivation. In its full shape, however, our theory provides means to follow the influence of finite temperatures, electron-electron scattering, and the reshaping of the exciton structure due to the presence of photons (such as forming the “flexible” exciton where many s-states are mixed) on the exciton nonlinearities. 

To provide possible explanation of experimentally-observed saturation constants, the dark exciton states are self-consistently introduced in the exciton-photon action via consideration of the electron and hole spins. 


The paper is organized as follows. Section <ref> is devoted to the full derivation of the exciton-photon action at finite temperature in the spinless case. Sec. <ref> addresses the limit T=0 and rigid 1s-exciton state, which allows to compare results with the existing works. In Sec. <ref> we present the corrections appearing in the renormalization of the Rabi splitting due to contributions of dark excitons to the UP and LP dispersions shifts when the spins of underlying electron and holes are self-consistently taken into account. Sec. <ref> concludes our studies. The main text is supplemented with Appendices <ref>–<ref> containing some cumbersome calculation details.



§ THE EXCITON-PHOTON ACTION
 


In the presented derivation, we consider the most relevant case n_ ex a^2 ≪1. The opposite limit can be attributed to the formation of the exciton insulator proposed by Keldysh and Kopaev <cit.>, which for polariton systems was studied in Ref. <cit.>. We start with the action for the electron-hole-photon system so far neglecting the spins of particles, in terms of the field operators of electrons Ψ_c(v) in the conduction (valence) bands with effective masses m_c(v) and the dispersions ε_c(v)( k) = ±(E_ g/2 + ħ^2  k^2/2m_c(v)), and the field operator Ψ_ ph of photons in a microcavity, with the effective mass m_ ph and the dispersion E_ ph( k) = E_ ph^0 + ħ^2  k^2/2 m_ ph. The origin of the energy is taken in the middle of the bandgap of width E_ g, and E_ ph^0 represents the cavity cutoff. The action reads:

    𝒮[Ψ_c,Ψ_v, Ψ_ ph] = ∫ d r∫_0^βdτ
    [[ Ψ_c(x) Ψ_v(x) ][ ∂_τ+ε_c(k̂) -μ_c                0;                0 ∂_τ+ε_v(k̂) -μ_v ][ Ψ_c(x); Ψ_v(x) ].  
    
            +  Ψ_ ph(x)(∂_τ + E_ ph(k̂) - μ_ ph) Ψ_ ph(x) 
    . + g_ R(Ψ_ ph(x)Ψ_c(x)Ψ_v(x) +Ψ_c(x)Ψ_v(x)Ψ_ ph(x) )] 
    
            + 1/2∑_i,j∫ d rd r^'∫_0^β dτ dτ^' V(x-x^') Ψ_i(x)Ψ_i(x) Ψ_j(x^') Ψ_j(x^'),

with i,j running over {c,v}, x = ( r, τ), τ is the imaginary time, β=ħ/k_ BT, k̂ is the momentum operator, and V(x-x^') ≡ V( r-r^')δ(τ-τ^'), where V( r-r^') represents the potential of interaction between charged particles. μ_c(v) is the chemical potential of the electrons in conduction (valence) band (which approximately equals zero for the ground state in undoped semiconductors), 
and μ_ ph is the chemical potential of photons. We consider separate chemical potentials of species since the total number of particles (electrons, holes, and photons) is not conserved. Finally, g_ R is the amplitude of light-matter coupling (electron-hole annihilation with photon creation and vice versa). In the functional integration (path integration) approach the grand canonical partition function is defined as:

    𝒵 = 𝒟[Ψ_c,v, Ψ_c,v]𝒟[Ψ_ ph, Ψ_ ph]e^-𝒮[Ψ_c,Ψ_v, Ψ_ ph] .

The exciton field Δ( r_1, r_2; τ) = Ψ_v( r_2,τ)Ψ_c( r_1,τ) can be introduced in the following way:

    𝒵 =    ∫𝒟[Ψ] e^-𝒮[Ψ]
    
        =    ∫𝒟[Ψ] 𝒟[Δ̅Δ] e^-𝒮[Ψ]δ(Δ - Ψ_vΨ_c) δ(Δ̅ - Ψ_cΨ_v) 
    
        =    ∫𝒟[Ψ] 𝒟[Δ̅Δ] 𝒟[ϕ̅,ϕ]e^-𝒮[Ψ]e^iϕ̅(Δ - Ψ_vΨ_c) + iϕ (Δ̅ - Ψ_cΨ_v),

where ϕ( r_1, r_2; τ) is an auxiliary field. Looking ahead, this field is unphysical and below will be integrated over.
Strictly speaking, other fields can be taken into account using the standard Hubbard-Stratonovich transformation: one can consider different channels of dual-field decoupling in the interacting part of the initial action. In particular, introducing the dual field Ψ_i(x)Ψ_i(y) (i=c,v), one can derive the renormalization of the exciton-exciton interaction matrix element due to electron scattering, which is done in Appendix <ref>. 

For convenience we perform the Fourier transform assuming the system homogeneous and finite, and account for finite temperatures:

    Ψ_i(x,τ)    =   √(T/S)∑_ k ∑_ω_nΨ_i( k, ω_n) e^i  k· x - iω_nτ,
    Δ(x, y)    =   √(T)/S∑_ k, k^'∑_Ω_nΔ( k,  k^', Ω_n) e^i k· x - i  k^'· y - iΩ_nτ,
    ϕ(x, y)    =   √(T)/S∑_ k, k^'∑_Ω_nϕ( k,  k^', Ω_n) e^i k· x - i  k^'· y - iΩ_nτ,

where ω_n=(2π n+π)/β and Ω_n=2π n/β are the fermionic and bosonic Mastubara frequencies, respectively, n∈ℤ, k, k^' are the quantized wavevectors, and S is the system area. The action takes the following form:


    𝒮[Ψ_c,Ψ_v,Ψ_ ph, ϕ, Δ]=∑_ k,ω[
            [ Ψ_c(k) Ψ_v(k) ][ -i ω + ε_c( k) - μ_c                    0;                    0 -i ω + ε_v( k) - μ_v ][ Ψ_c(k); Ψ_v(k) ]] 
    
        + i√(T)∑_ k_1,ω_1∑_ k_2,ω_2[
            [ Ψ_c(k_1) Ψ_v(k_1) ][                      0  ϕ( k_1, k_2, ω_1-ω_2); ϕ̅( k_2, k_1, ω_2-ω_1)                      0 ][ Ψ_c(k_2); Ψ_v(k_2) ]] 
    
        - i ∑_ k_1, k_2,Ω[ ϕ̅( k_1, k_2, Ω)Δ( k_1, k_2, Ω) + Δ̅( k_1, k_2, Ω) ϕ( k_1, k_2, Ω) ]
        
    
        - ∑_ k_1… k_4, Ω V( k_1- k_2) Δ̅( k_1, k_4,Ω) Δ( k_2, k_3,Ω) δ( k_1 + k_3,  k_2 + k_4) + ∑_ k, ΩΨ_ ph(k)(- i Ω + E_ ph( k) - μ_ ph)Ψ_ ph(k) 
    
        + g_ R/√(S)∑_ k_1, k_2,Ω[Ψ_ ph( k_1 - k_2,Ω) Δ( k_1, k_2,Ω) + Δ̅( k_1, k_2,Ω)Ψ_ ph( k_1 - k_2,Ω)],


where all sums over ω,Ω involve summation over integer indices n (subscripts n are omitted for clarity), and the four-vector notation k ≡ ( k, ω_n) or ( k,Ω_n) was introduced for both fermionic and bosonic fields.

It is then straightforward to perform the integration over the fermionic fields Ψ_c,v, which results in the (- Trln𝒢^-1) contribution to the action, where 𝒢^-1 = 𝒢^-1_0+ δ𝒢^-1 with 

    𝒢^-1_0=[ -i ω_1 +ε_c( k_1) -μ_c                      0;                      0 -i ω_1 +ε_v( k_1) -μ_v ]δ_k_1k_2

denoting the bare electron Green's function and 

    δ𝒢^-1= i√(T)[                             0   ϕ( k_1, k_2, 
               Ω); ϕ̅( k_2, k_1, 
               -Ω)                             0 ]

containing the introduced auxiliary field ϕ (here the bosonic Matsubara frequency Ω≡Ω_n originates from the difference of two fermionic frequencies ω_1n - ω_2n). It can be easily calculated under the assumption δ𝒢^-1≪𝒢^-1_0 (considering δ𝒢^-1 as a perturbation, which corresponds to our applicability condition n_ ex a^2 ≪1):

    Trln𝒢^-1 =  Trln𝒢_0^-1 +  Trln(1+𝒢_0δ𝒢^-1) =
    =   Trln𝒢_0^-1 +  Tr∑_n(-1)^n-1/n(𝒢_0δ𝒢^-1)^n.

In the expansion series one should keep the terms up to the 4th power to derive the exciton-exciton interaction. 

Since 𝒢_0 is diagonal while δ𝒢^-1 is antidiagonal, the first and third terms in the logarithm expansion are absent, hence in the lowest order in the auxiliary field ϕ one gets:

    1/2 Tr(𝒢_0δ𝒢^-1)^2 =     -T ∑_ k_1,  k_2, Ωϕ̅( k_1, k_2, Ω)ϕ( k_1, k_2, Ω) 
       ×∑_ω^'𝒢_0( k_1, ω^'+Ω)_11𝒢_0( k_2, ω')_22.

The last summation is performed over the fermionic Matsubara frequency ω^' of the loop diagram, and for finite temperatures it leads to the exciton propagator

    T ∑_ω^'𝒢_0( k_1, ω^' + Ω)_11𝒢_0( k_2, ω^')_22 = 
    
      =n_v( k_2)-n_c( k_1)/ε_c( k_1) -ε_v( k_2) -μ_c+μ_v- iΩ≡𝒜_ k_1 k_2^  Ω,

where n_i( k) = 1/(e^β (ε_i( k) -μ_i)+1) is the Fermi distribution function in the i-band  (the details are provided in Appendix <ref>). Here we note that for an undoped semiconductor n_c = n_v and thus one can approximate μ_c∼ -μ_v + Tln(m_v/m_c). From the point of view of equilibrium chemical reactions, the exciton chemical potential is defined as μ_ ex = μ_e + μ_h = μ_c - μ_v (where μ_e,h denote the chemical potentials of the electrons and holes), and in thermodynamic equilibrium between the excitons and photons 
it should equal μ_ ex = μ_ ph≡μ. Then the above expression reduces to

    𝒜_ k_1 k_2^  Ω =  n_v( k_2)-n_c( k_1)/ε_c( k_1) - ε_v( k_2) -μ-i Ω.


In the same manner we calculate the |ϕ|^4-term:

    1/4 Tr(𝒢_0δ𝒢^-1)^4 = 2/4(i√(T))^4∑_ k_1… k_4∑_Ω_1…Ω_3ϕ( k_1, k_2, Ω_1)
    ×ϕ̅( k_3,  k_2, Ω_2)ϕ( k_3,  k_4, Ω_3)ϕ̅( k_1,  k_4, Ω_1 +Ω_3 -Ω_2)
        
    ×∑_ω^'𝒢_0( k_1, ω^')_11𝒢_0( k_2, ω^'-Ω_1)_22𝒢_0( k_3, ω^'-Ω_1 +Ω_2)_11
    ×𝒢_0( k_4, ω^'- Ω_1+ Ω_2- Ω_3)_22,

where the factor of 2 arises from the number of loops in the corresponding diagram series (see Appendix <ref> and Fig. <ref> therein). 
The result of summation over the frequency of the fourth-order loop diagram in Eq. (<ref>) is given in Appendix <ref>, while here we only introduce the loop notation:


    ℒ_ k_1 k_2 k_3 k_4^Ω_1Ω_2Ω_3 = T ∑_ω^'𝒢_0( k_1, ω^')_11𝒢_0( k_2, ω^'-Ω_1)_22
    ×𝒢_0( k_3, ω^'-Ω_1 +Ω_2)_11𝒢_0( k_4, ω^' -Ω_1 +Ω_2 -Ω_3)_22

The effective action takes the form:


    𝒮[Ψ_ ph,ϕ, Δ]= - ∑_ k_1,  k_2, Ω𝒜_ k_1 k_2^  Ω ϕ̅( k_1, k_2, Ω)ϕ( k_1, k_2, Ω) + ∑_ k, ΩΨ_ ph(k)(- i Ω + E_ ph( k) -μ)Ψ_ ph(k)
         
    
         -i ∑_ k_1, k_2,Ω[ϕ̅( k_1, k_2, Ω)Δ( k_1, k_2, Ω) + Δ̅( k_1, k_2, Ω) ϕ( k_1, k_2, Ω) ]
         - ∑_ k_1… k_4, Ω V( k_1 -  k_2) Δ̅( k_1, k_4,Ω) Δ( k_2, k_3,Ω) δ( k_1+ k_3,  k_2+ k_4)
         
    
         + g_ R/√(S)∑_ k_1, k_2,Ω[Ψ_ ph( k_1- k_2,Ω) Δ( k_1, k_2,Ω) + Δ̅( k_1, k_2,Ω)Ψ_ ph( k_1- k_2,Ω)]
         
    
         +T/2∑_{ k_i, Ω_i}ℒ_ k_1 k_2 k_3 k_4^Ω_1Ω_2Ω_3 ϕ( k_1, k_2, Ω_1)ϕ̅( k_3,  k_2, Ω_2)ϕ( k_3,  k_4, Ω_3)ϕ̅( k_1,  k_4, Ω_1 + Ω_3 -Ω_2).



Now, the auxiliary field ϕ needs to be excluded. Within the saddle-point approximation in the lowest order of the Trln𝒢^-1 expansion we have:

    δ𝒮/δϕ̅= - 𝒜_ k_1 k_2^  Ω ϕ( k_1, k_2,Ω) - i Δ( k_1, k_2,Ω) =0,

and the relation between the exciton and the auxiliary fields can be found:

    
    ϕ( k_1, k_2,Ω) = Δ( k_1, k_2,Ω)/i𝒜_ k_1 k_2^  Ω, 
    ϕ̅( k_1, k_2,Ω) =Δ̅( k_1, k_2,Ω)/i𝒜_ k_1 k_2^  Ω,


which allows one to rewrite the effective action in terms of the exciton and photon fields, and obtain the effective exciton-photon action. Despite the result being quite similar with the previously reported <cit.>, in our case the exciton field depends on two momenta. We can interpret it as the exciton field Δ( k_1, k_2, Ω) corresponding to the relative motion of electrons belonging to different bands, with the momentum p = (m_c k_1+m_v k_2)/(m_c+ m_v), and motion of the exciton as a whole, with the momentum k= k_1- k_2.



 §.§ Bethe-Salpeter equation for excitons coupled to photons at finite temperature




It is useful to rewrite the obtained expressions in terms of the relative and total momenta:

    𝒮[Ψ_ ph,Δ]= ∑_ p,  k, ΩΔ̅( p, k, Ω)Δ( p, k, Ω)/𝒜_ p k^ Ω
    
          + ∑_ k, ΩΨ_ ph( k,Ω)(- i Ω + E_ ph( k) -μ)Ψ_ ph( k,Ω) 
    
          - ∑_ k, p,  q, Ω V( p - q) Δ̅( p, k,Ω) Δ( q, k,Ω)  
    
          +g_ R/√(S)∑_ p, k,Ω[Ψ_ ph( k,Ω) Δ( p, k,Ω) +Δ̅( p, k,Ω)Ψ_ ph( k,Ω)]
    
          + T/2∑_ p,{ l_i, Ω_i}ℒ̃_ p  l_1 l_2 l_3^Ω_1 Ω_2Ω_3 Δ( p, l_1, Ω_1)Δ̅( p + l_2 - l_12,  l_2, Ω_2) 
    ×Δ( p + l_2-  l_1 - l_32,  l_3, Ω_3)Δ̅( p +  l_2 - l_32,  l_1 -  l_2+  l_3, Ω_1 + Ω_3 - Ω_2).

Here and below, tildes denote the loops connecting exciton fields Δ (see Appendix <ref>). When performing the saddle-point approximation to define the shape of the exciton field in presence of photons, the last term can be considered as perturbation since it was derived in the higher order in the expansion series, so that 
it can be temporarily neglected:


    δ𝒮/δΔ̅   =   Δ( p, k, Ω)/𝒜_ p k^ Ω + g_ R/√(S)Ψ_ ph( k, Ω) 
             - ∑_ qV( p - q) Δ(  q, k,Ω)=0, 
    δ𝒮/δΨ_ ph   =    (- i Ω + E_ ph( k) -μ)Ψ_ ph( k,Ω) 
                 + ∑_ pg_ R/√(S)Δ( p, k, Ω)=0.


Determining the photon field via the exciton field in (<ref>) and substituting it into (<ref>), we obtain the generalized Bethe-Salpeter equation modified by the presence of the photon field (similar to the semiconductor Bloch equation <cit.> albeit with the explicit dependence on k in the photon dispersion):

    Δ( p, k, Ω)/𝒜_ p k^ Ω - ∑_ q[V( p - q) . 
    
        + . (g_ R/√(S))^21/E_ ph( k) -μ- i Ω]Δ(  q, k,Ω)=0.


This equation is one of the main results of this study. Finite temperatures are contained in the summations and Fermi distributions appearing in 𝒜_ p k^ Ω. As anticipated, the photon field alters the exciton interaction, since the electrons and holes in such a system can be bound not only by Coulomb-like interactions, but also by light. This phenomenon is known in the literature as the flexible exciton limit <cit.> in which the polariton Rabi splitting is not negligible compared to the exciton binding energy: ħΩ_ R∼ E_ b.
It should be noted that similar equations were recently derived in Ref. <cit.> for the zero-temperature case using the operator formalism. Ref. <cit.> proposes the solution of (<ref>) for T=0 in the case of one-mode fields. To compare, we apply this method to the equation (<ref>) at zero temperature in Appendix <ref>.

Eq. (<ref>) can be solved using the ansatz 
    Δ( p,  k,Ω) = ∑_νχ^(ν)_ k, Ω( p)C^(ν)( k, Ω),

which is a decomposition of the exciton field over the basis of (hydrogen-like) wavefunctions χ^(ν)_ k,Ω( p), with ν being the number of excitonic ν s–state and C^(ν)( k, Ω) the corresponding exciton field depending only on the total momentum. Introducing the interaction matrix element g( l_1, l_2, l_3, Ω_1,Ω_2,Ω_3,ν_1,ν_2,ν_3,ν_4)≡ g( l_i, Ω_i, ν_i) as

    g( l_i, Ω_i, ν_i) =∑_ pℒ̃_ p  l_1 l_2 l_3^Ω_1 Ω_2Ω_3 χ^(ν_1)_ l_1, Ω_1( p)χ̅^(ν_2)_ l_2, Ω_2( p + l_2 - l_12)
        
        ×χ^(ν_3)_ l_3, Ω_3( p + l_2-  l_1 - l_32)χ̅^(ν_4)_ l_1-  l_2+  l_3, Ω_1-Ω_2+ Ω_3( p +  l_2 - l_32),

we arrive at the final expression for effective exciton-photon action:

    𝒮[Ψ_ ph, C]= ∑_ k, Ω∑_νC^(ν)( k, Ω)E^(ν)_ ex( k, Ω) C^(ν)( k, Ω)
         
    + ∑_ k, ΩΨ_ ph( k,Ω)(- i Ω + E_ ph( k) -μ)Ψ_ ph( k,Ω) 
    
         + g_ R∑_ k, Ω∑_ν[χ̅^(ν)_ k, Ω( r= 0) C^(ν)( k, Ω)Ψ_ ph( k, Ω) 
         + 
         c.c.]
    
        + T/2∑_{ l_i, Ω_i, ν_i} g( l_i, Ω_i, ν_i)C^(ν_1)( l_1, Ω_1)C^(ν_2)(  l_2, Ω_2)
    × C^(ν_3)( l_3, Ω_3)C^(ν_4)(  l_1-  l_2+  l_3, Ω_1 + Ω_3 -Ω_2).

Here it is important to note that all the expressions are derived in the most general form (for T 0, arbitrary interactions and for excitons and photons with k 0). Summation over ν allows to account for the flexible exciton case when the coupling to light results in the mixing of states in the exciton s-series. Appendix <ref> is devoted to simplifications of these general formulae for the case of zero temperature.



 §.§ Existing limits



Standard treatment of exciton-polariton systems often neglects 
the second term in the square brackets in Eq. (<ref>) or, analogously, neglects the terms ∼ g_ R compared to the exciton binding energy E_ b in Eqs. (<ref>,<ref>). Such considerations correspond to the rigid exciton limit when E_ b≫ħΩ_ R:

    1/𝒜_ p k^ ΩΔ( p, k, Ω) - ∑_ qV( p - q) Δ(  q, k,Ω)=0.

From this point of view, considering only the 1s-state at T=0  and separating variables Δ( p,  k,Ω) = χ( p)C( k,Ω) (where χ( p)≡χ^(1)_ k,Ω( p) is independent of k and Ω) brings Eq. (<ref>) to the standard Wannier equation. In this limit, one would obtain the exciton wavefunction and the renormalization of the exciton-photon conversion term previously reported in several works <cit.>:

    g_ R/√(S)∑_ qχ( q) = g_ Rχ( r = 0) = ħΩ_ R/2,

where ħΩ_ R/2 is the experimentally-observed polariton Rabi-splitting. In this case, the saddle-point equations (<ref>, <ref>) become the standard Hopfield equations while the exciton-photon action takes the well-known form <cit.>. Here we emphasize that this approach to treat the action is, in principle, equivalent to the assumption that first an electron and a hole are bound by electrostatic interaction forming an exciton, and then excitons are coupled to photons, which is a common way to describe strong coupling between excitons and photons. Strictly speaking, in the case of finite temperatures variables cannot be divided, as can be anticipated from the shape of 𝒜_ p k^ Ω in Eq. (<ref>) with ( k_1, k_2)→( p, k).

On the other hand, in the opposite limit of dominating photon-electron-hole coupling, Eq. (<ref>) for the exciton field Δ is similar to (<ref>) but with the electron-hole interactions induced only by light <cit.>, which means one cannot parametrize the field Δ as χ( p)C( k, Ω) with the hydrogen-like wavefunction χ( p). 



§ INTERACTION CONSTANTS
 


To describe the nonlinearities in terms of interaction constants including the so-called saturation (exciton-assisted photon-exciton coupling) that arises in our treatment from the fourth-order expansion term, one needs not only to set T=0 but also to assume the rigid 1s-exciton limit. While for the flexible exciton where all ν s states of the hydrogenic basis get mixed, the saturation contribution is also derived (see Appendix <ref>), it cannot be presented as a closed-form explicit expression. The case of rigid excitons, in which we can analyze the results of recent experiments <cit.> and compare to existing theoretical results <cit.>, is considered in Appendix <ref>.

For the ground-state exciton wavefunction which in momentum space reads χ( p) = 2 a√(2π)/(1 + p^2 a^2)^3/2, the expression for the saturation interaction constant is defined as g_ sat = (8 π/7) a^2ħΩ_ R/2, where a in the hydrogenic description is the 2D exciton Bohr radius, while for the Rytova-Keldysh 1s-state should be considered a variational parameter <cit.>. This is a well-established result obtained in many different approaches <cit.>. Clearly, the exciton-assisted exciton-photon coupling ∼ g_ satn_ ex describes the conversion between photons and excitons in the higher order of the small parameter n_ exa^2 compared to the terms ∼(Ψ_ ph C +  c.c.). 

The corresponding reduction of the polariton Rabi splitting can be defined from the shifts of the UP and LP dispersions with the increase of the exciton density. Performing the saddle-point approximation for excitons and photons in the assumption of uniform exciton density n_ ex = const (see Appendix <ref>), we obtain 






the splitting between the dispersions minima renormalized by interactions 

    ħΩ_R(n_ ex)   =   ħΩ_ R√(1  - 4 g_ satn_ ex/ħΩ_ R + 3 (g_ satn_ ex)^2 + g_ exn_ ex/(ħΩ_R)^2)
       ≈   ħΩ_ R - 2g_ satn_ ex.

Since g_ satn_ ex∼ħΩ_ R n_ exa^2 and the exciton densities considered here are much smaller than the Mott density ∼ a^-2, the correction to the splitting is supposed to be much smaller than its value without the renormalization.

Futhermore, we show (see Appendices <ref> and <ref>) that in the rigid 1s-exciton limit the pure exciton-exciton contribution to interaction for Wannier-Mott excitons should be larger than that from saturation, since 
the Rabi-coupling term in the saddle-point Eqs. (<ref>,<ref>) in this limit is treated as perturbation. As mentioned in the Introduction, the exciton interaction constant in this limit at T=0 for excitons in QWs (i.e. with Coulomb interaction potential) is equal g_ ex= (1 - 315π^2/4096)4πħ^2/M=6.06 E_ b a^2 <cit.>. 

For TMD-based materials (with Rytova-Keldysh interaction potential) at T=0 we derive from Eq. (<ref>) g_ ex = (16/π) e^2 a f(r_0/a), where the overlap f(r_0/a) of the 1s–exciton wavefunction and Fourier image of the interaction potential depends only on the ratio between the screening length r_0 and the Bohr radius a. Importantly, due to the reduced screening and hence small spread of the wavefunction, this overlap is small and decreasing with the growth of r_0. As a consequence, the exciton interaction constant g_ ex for 2D materials with the Rytova-Keldysh interaction between the charges is smaller than that for the 2D Coulomb potential.

It is important to underline that, as can be seen from Eqs. (<ref>) and (<ref>) even without the substitution of χ( p) in the hydrogenic shape, the constants g_ ex and g_ sat are defined solely by the properties of the material and cannot be used as tunable fitting parameters, unless one assumes that the rigid exciton limit is violated and thus the theoretical expressions obtained in this limit are not valid. In this light, the conclusions drawn in Ref. <cit.> from direct measurements of the Rabi splitting at different exciton densities in monolayer WS_2 that saturation in their system is huge, requires reconsideration. Furthermore, the recent experiment <cit.> performed for monolayer MoSe_2 revealed from fitting g_ sat≫ g_ ex (in <cit.>, g_ ex is totally disregarded). However, since the excitons in TMD materials are robust (E_ b≫ħΩ_ R), the saturation process should be a correction rather than the leading term. We note that for the case of GaAs QWs <cit.> where E_ b≳ħΩ_ R, the exciton rigidity can be undermined by the presence of photons and g_ sat could be of the same order with g_ ex, but the effect of flexibility (mixing of the excitonic s-series) in TMDs is highly unlikely.


In Table <ref>, we summarize the existing experimental estimates for g_ ex and g_ sat obtained from fitting of the polariton branches in TMD monolayers, and compare them to each other and to the theoretical values obtained here. One sees that only the experimental values obtained at low temperature (5 K) in Ref <cit.> have the same order as the theoretically-predicted constants. Hence, noting that straightforward introduction of the interaction constants works only at T→ 0 when one can factorize the exciton field (since χ_ k( p) stops being dependent on Ω), one needs to rely on the full temperature-dependent treatment [see Eq. (<ref>)]. More strikingly, the nonlinearities experimentally reported for the same materials in different works (even at the same T) differ by orders of magnitude. We conclude that while the rigid-exciton limit for TMD-materials is expected to be valid with a great accuracy, there are other mechanisms at play strongly altering the polariton nonlinearities. 




§ RENORMALIZATION OF THE RABI SPLITTING DUE TO DARK EXCITONS
 


In this section, we show that one of such possible mechanisms to explain this inconsistency can be addressed by taking into account the spins of particles and considering the influence of dark states on the saturation process. For simplicity, we restrict ourself only to the case of zero temperature. The details of derivations, which in spirit are analogous to those described in Sec. <ref>, are presented in Appendix <ref>. In this case, four excitonic fields are introduced, two of which correspond to bright Δ_± 1 and the other two to dark excitons Δ_ d,d̅ (here d, d̅ generically denote the two spin projections of excitons that are not coupling to light, i.e. ±2 for GaAs and 0 for TMD excitons in two different valleys). Moreover, in Appendix <ref> we take into account possible spin-bands splitting, i.e. different dispersion curvatures for electrons with different spin projections. The resulting saddle-point equations for the four excitonic fields are differ from each other due to the deviation in the exciton dispersions and, consequently, in the 1s–exciton wavefunctions. As a result, the exciton interaction and saturation constants for Δ_±1 and Δ_ d,d̅ generally differ as well. 

Here, to show the simplest case, we assume that photons of only one polarization (either +1 or -1) are present in the system and that the effective electron masses m_c(v) corresponding to different spin-bands are equal. Then one obtains from the general expression the splitting between the UP and LP branches at k=0:

    ħΩ_R(n_ ex, n_ d) = ħΩ_ R[1  -  g_ sat(4 n_ ex + 2 n_ d)/ħΩ_R.
     + .g^2_ sat(n_ ex+ n_ d)(3n_ ex+ n_ d) + g^2_ ex(n_ ex+ n_ d)^2/(ħΩ_ R)^2]^1/2,

where n_ d = C_d̅ C_d̅+ C_ d C_ d is the dark exciton density, C_ d(d̅) is the dark-exciton field depending only on the total exciton momentum. Strictly speaking, in Eq. (<ref>) we took into account only the dominant terms (for details see Appendix <ref>). Surprisingly, n_ ex and n_ d enter this expression additively with the coefficients of the same order, which underlines that the dark exciton population cannot be disregarded when addressing the exciton saturation. The presence of dark exciton states, as clearly seen, leads to further quench of the Rabi splitting with the growth of both densities, which could be the possible explanation of the abovementioned experimental observations. It is worth noting that a similar model including the reservoir contribution was  phenomenologically proposed in Ref. <cit.>, while the reservoir was assumed to result from biexciton states. 
Our model, however, does not take the exciton bound states into account, thus overlooking the interactions of excitons with antiparallel spins g_ ex^↑↓ |C_+1|^2|C_-1|^2 which can contribute to polariton interactions <cit.>. In the existing literature <cit.> the corresponding interaction constant g_ ex^↑↓ is demonstrated to be negative and, in general, comparable with g_ ex. Nevertheless, on the qualitative level, these attractive interactions are expected to cause the decrease of blue- and redshifts of the polariton dispersions. 



§ CONCLUSIONS
 


To summarize, we developed the approach to perturbatively describe the electron-hole-photon system in presence of strong coupling, based on the equilibrium path integral technique. This approach, in general, allows to track different channels of dual-field pairing, and self-consistently takes into account finite temperatures and the full exciton ν s-series. It can be applied to any material system (such as conventional quantum-well microcavities or those based on TMD monolayers and bilayers). The difference in the charges interaction potentials will enter via the shape of the exciton wavefunctions that define the final quantitative results. As a logical extension we also provided (see Appendix <ref>) the treatment of particle spins, which revealed the contributions to exciton nonlinearities due to bright-bright, bright-dark, dark-dark exciton interactions and spin-flip processes. 

In the simplified case of rigid excitons, our analysis provides important results relevant to understanding the giant nonlinearities observed in TMD-based microcavities. First, we show that finite temperatures result in a sizeable change of the interaction constants and the usual zero-T expressions cannot be used for the estimates in room-temperature studies. Second, we reveal that the exciton-photon conversion is affected by the presence of the reservoir particles (such as dark excitons). We rigorously derive the terms involving the dark exciton fields in effective action, and the dark-exciton density contributions to the Rabi splitting renormalization. These results call for a detailed  comparative experimental studies of the exciton saturation in TMD structures, to be performed at different temperatures, which would allow to elucidate the relative influence of the two factors on the observed nonlinearity.

In a more sophisticated case of flexible excitons, which is more relevant to quantum-well polaritons or for situations with very strong coupling, an important conclusion is that the concept of interaction constants describing the two sources of nonlinearity in the system starts to fail. As several exciton s-states start to contribute to the exciton field, the relative motion of the electron and hole within the exciton cannot be separated from the motion of the exciton as a whole. In this case the nonlinearities arising from the fourth-order term with respect to the exciton field are not factorized to the integral over the exciton fields and wavefunctions overlap integrals (which yield the above mentioned interaction and saturation constants). Furthermore, the Bethe-Salpeter equation in this case contains not only the temperature dependence, but also the contribution of the momentum-dependent photon energy dispersion to the formation of the exciton.

We note that our study did not consider additional contributions to the photon-exciton action such as the formation of multi-exciton bound states. These scenarios were addressed in recent works <cit.> where it was shown that biexciton-polariton formation is possible, and that biexciton transition therefore can also alter nonlinear optical response of TMDs.



§ ACKNOWLEDGEMENTS



The authors are thankful to Vanik Shahnazaryan for useful discussions. The partial financial support from the NRNU MEPhI Priority 2030 Program is acknowledged. The work of A.G. is funded by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”, under the Grant No. 22–1–5–30–1.




§ ADDITIONAL FIELDS (OTHER CHANNELS OF PAIRING)

In the beginning of the main text we introduce the exciton field only. Strictly speaking, 
such consideration overlooks the exciton interaction mediated by electron scattering processes. In this Appendix, we address this problem and show how taking into account the 
off-diagonal density field Ψ_i( r, τ)Ψ_i( r^', τ) (i=c,v) changes the exciton interaction term. For clarity, in the derivations of this section we will omit the “photon” part of Eq. (<ref>) which does not play a role here, being not coupled with the density fields.
Turning to the initial form of the action (<ref>), we introduce not only the previously discussed fields Δ and ϕ, but also other density-like fields which arise when considering other pairing channels, such as the off-diagonal Φ( r,  r^', τ) = Ψ_i( r, τ)Ψ_i( r^', τ) and diagonal ξ_i( r, τ) = Ψ_i( r, τ)Ψ_i( r, τ) density fields, using the Hubbard-Stratonovich transformation <cit.>. 

First, we define the momentum-frequency representation for the considered fields:

    Φ_i(x, y) = √(T)/S∑_ k, k^'∑_  Ω_nΦ_i( k,  k^', Ω_n) e^i k· x - i  k^'· y - iΩ_nτ
    ξ_i(x) = √(T/S)∑_ k, k^'∑_  Ω_nξ_i( k,  k^', Ω_n) e^i k· x - i k^'· x- iΩ_nτ,

where k are discrete momenta and Ω_n  = 2π n T ≡Ω are bosonic Matsubara frequencies. As one can see, the Fourier image of the diagonal density field depends on one momentum (k- k^') only, which leads, as we will see further, to qualitatively different contributions to the exciton interaction.
In the frequency-momentum representation the action takes the form:

    𝒮[Ψ_c,Ψ_v, ϕ, Δ, Φ_c, Φ_v, ξ_c, ξ_v]=∑_ k,ω[
           [ Ψ_ c(k) Ψ_ v(k) ][ -i ω + ε_c( k) - μ_c                    0;                    0 -i ω + ε_v( k) - μ_v ][ Ψ_ c(k); Ψ_ v(k) ]]


    -i ∑_ k_1, k_2
    Ω[ ϕ̅( k_1, k_2, Ω)Δ( k_1, k_2, Ω) + Δ̅( k_1, k_2, Ω) ϕ( k_1, k_2, Ω) ] - ∑_ k_1… k_4
    Ω V( k_1 -  k_2) Δ̅( k_1, k_4,Ω) Δ( k_2, k_3,Ω) δ( k_1+ k_3,  k_2+ k_4)
        
    
        + 1/2∑_ k, Ω
        V^-1( k)ξ( k, Ω)ξ(- k, -Ω) + 1/2∑_ k_1… k_4
    Ω∑_  i=c,v V^-1( k_1 -  k_2) Φ_i( k_1, k_4,Ω) Φ_i( k_3, k_2,Ω) δ( k_1+ k_3,  k_2+ k_4) 
    
        + √(T)∑_ k_1,ω_1∑_ k_2,ω_2[
            [ Ψ_ c(k_1) Ψ_ v(k_1) ][ Φ_c(k_2, k_1) + iξ(k_2- k_1)    i ϕ( k_1, k_2, ω_1 - ω_2);   i ϕ̅( k_2, k_1, ω_2 - ω_1) Φ_v(k_2,k_1) +  iξ(k_2- k_1) ][ Ψ_ c(k_2); Ψ_ v(k_2) ]] + …,

where, compared to Eq. (<ref>) the terms containing Ψ_ ph are omitted (…) and the notation ξ = ξ_c + ξ_v is introduced.
Although the expression is quite cumbersome, all the contributions containing Φ and ξ will be integrated out, and the final expression remains elegant.
After performing the path integration over the fermionic fields, the second order of the Tr ln𝒢^-1 expansion series yields the relation between ϕ and Δ, as in the main text,

    ϕ( k_1, k_2,Ω) = -i/𝒜_ k_1 k_2^  ΩΔ( k_1, k_2,Ω),

and renormalizes the Coulomb interaction for the 
off-diagonal and diagonal densities

    V^-1( k_1 - k_2)    →    V^-1( k_1 - k_2) + T ∑_ω^'𝒢_0( k_1, ω^'+ Ω)_11𝒢_0( k_2, ω^')_11  for Φ_c, 
    
        V^-1( k_1 - k_2)    →    V^-1( k_1 - k_2) + T ∑_ω^'𝒢_0( k_1, ω^'+ Ω)_22𝒢_0( k_2, ω^')_22  for Φ_v, 
    
        V^-1( k_1 - k_2)    →    V^-1( k_1 - k_2) - T ∑_ω^'𝒢_0( k_1, ω^'+ Ω)_11𝒢_0( k_2, ω^')_11+ T∑_ω^'𝒢_0( k_1, ω^'+ Ω)_22𝒢_0( k_2, ω^')_22  for ξ,

where ω^' = (2n+1)π T is the fermionic Matsubara frequency. 

The renormalization (<ref>) is assumed to be small, since we assume the small-parameter expansion. We note that since the variation of the Green's function δ𝒢^-1 contains both diagonal and off-diagonal elements, the first- and third-order terms appear in the expansion series. The first-order term is zero under the assumption that the system is electrically neutral. 
The third-order term has the form: 

    1/3 Tr(𝒢_0δ𝒢^-1)^3 = 1/3∑_k_1, k_2, k_3𝒢_0(k_1)_ααδ𝒢^-1(k_1,k_2)_αβ𝒢_0(k_2)_ββδ𝒢^-1(k_2, k_3)_βγ𝒢_0(k_3)_γγδ𝒢^-1(k_3, k_1)_γα
    
        = √(T)∑_ k_1, k_2, k_3∑_Ω_1, Ω_2[Π̃_121]_ k_1 k_2 k_3^ Ω_1 Ω_2 Δ( k_1, k_2, Ω_1) Δ̅( k_3,  k_2, Ω_2) [Φ_c( k_1,  k_3, Ω_1 - Ω_2)+ iξ( k_1-  k_3, Ω_1 - Ω_2)] 
    
        + √(T)∑_ k_1, k_2, k_3∑_Ω_1, Ω_2[Π̃_122]_ k_1 k_2 k_3^ Ω_1 Ω_2 Δ( k_1, k_2, Ω_1) Δ̅( k_1,  k_3, Ω_2) [Φ_v( k_3,  k_2, Ω_1 - Ω_2) + i ξ( k_3-  k_2, Ω_1 - Ω_2)],

where α,β,γ run over the values 1,2, but we choose only those that contain the fields ϕ (or Δ), hence the notation for the third-order loops Π_121 and Π_122 is introduced (see Appendix <ref> for definitions and details of calculation). Tilde, as before, denotes the correspondence of a loop to the exciton field Δ (i.e. after the exclusion of the auxiliary field ϕ).

Integrating over the fields Φ_c,v we obtain the correction to the exciton interaction due to the electron-electron scattering (i.e. screening):

    𝒮^'_ scr = - T/2∑_ q, k_1… k_4∑_Ω_1,Ω_2 
    Ω_1^',Ω_2^' V( q) [Π̃_121]_ k_1  k_3 k_2^ Ω_1Ω_2[Π̃_121]_ k_2 -  q,  k_4,  k_1 -  q^  Ω_1^'  Ω_2^' Δ( k_1, k_3, Ω_1)Δ̅( k_2, k_3, Ω_2)Δ( k_2 -  q,  k_4, Ω_1^') 
    ×Δ̅( k_1- q,  k_4, Ω_2^')δ_Ω_1 + Ω_1^' ,Ω_2 + Ω_2^'
        -T/2∑_ q, k_1… k_4∑_Ω_1, Ω_2 
    Ω_1^', Ω_2^' V( q) [Π̃_122]_ k_3  k_2 k_1^ Ω_1 Ω_2[Π̃_122]_ k_4,  k_1 -  q,  k_2 -  q^  Ω_1^'  Ω_2^'Δ( k_3, k_2, Ω_1)
    
       ×Δ̅( k_3, k_1, Ω_2)Δ( k_4,  k_1 -  q, Ω_1^')Δ̅( k_4,  k_2 -  q,Ω_2^') δ_Ω_1 + Ω_1^' ,Ω_2 + Ω_2^'.

Eq. (<ref>) has the same structure as the fourth-order term in the Trln𝒢^-1 expansion [see Eq. (<ref>)] and therefore is a correction to the |Δ|^4-term.

In a similar fashion, integrating out the diagonal density fields ξ_c,v leads to the |Δ|^4-term describing the repulsive exciton-exciton interaction 

    𝒮^'_ vdW = T/2∑_ q, k_1… k_4∑_Ω_1, Ω_2 
    Ω_1^',Ω_2^' V( q) [Π̃_121]_ k_1,  k_2,  k_1 -  q^  Ω_1  Ω_2[Π̃_121]_ k_3,  k_4,  k_3 +  q^  Ω_1^'  Ω_2^' Δ( k_1, k_2, Ω_1) Δ̅( k_1 -  q, k_2 , Ω_2)Δ( k_3, k_4, Ω_1^') 
    ×Δ̅( k_3 +  q, k_4, Ω_2^')δ_Ω_1 + Ω_1^', Ω_2 + Ω_2^'
       +T/2∑_ q,  k_1… k_4∑_Ω_1, Ω_2 
    Ω_1^',Ω_2^' V( q) [Π̃_122]_ k_1,  k_2,  k_2- q^  Ω_1   Ω_2[Π̃_122]_ k_3,  k_4,  k_4 +  q^  Ω_1^'  Ω_2^'Δ( k_1, k_2, Ω_1)
       
    
       ×Δ̅( k_1, k_2 -  q, Ω_2)Δ( k_3, k_4, Ω_1^') Δ̅( k_3, k_4  +  q, Ω_2^') δ_Ω_1 + Ω_1^',Ω_2 + Ω_2^'

which can be interpreted as a correction to the Coulomb-like interaction (-VΔ̅Δ), since 
it contains the direct Coulomb interaction of excitons as a whole. 
Therefore Eq. (<ref>) may be treated as a correction to the generalized Bethe-Salpeter equation due to the van-der-Waals interaction

    δ S/δΔ̅= 1/𝒜_ p k^ Ω Δ( p, k, Ω) - ∑_ q[V( p - q) + (g_ R/√(S))^21/E_ ph( k) -μ- i Ω]Δ(  q, k,Ω) + δ (𝒮^'_ vdW)/δΔ̅=0


which changes Eq. (<ref>) just slightly and can be neglected. We conclude that the 
consideration of additional pairing channels in the initial action (<ref>), which is performed here generally, leads to corrections due to the screening of the exciton-exciton interaction only.



§ CALCULATION OF THE LOOP DIAGRAMS


To derive the exciton propagator [the second-order loop shown in Fig. <ref>(a)], we use the standard resummation procedure, i.e. the Sommerfeld-Watson transformation:

    𝒜_ k_1 k_2^  Ω = T ∑_ω^' = (2n+1)π T𝒢_0( k_1, ω^'+ Ω)_11𝒢_0( k_2, ω^')_22 = T ∑_ω^'1/[-i(ω^'+ Ω) + ε_c( k_1) - μ_c][-iω^'+ ε_v( k_2)- μ_v ]
    
       =T1/2π i1/T∮ dz  1/e^β z+11/[-z -i Ω + ε_c( k_1)- μ_c][-z + ε_v( k_2)- μ_v]
    
       =1/ε_c( k_1) - ε_v( k_2) -μ_c + μ_v- i Ω(1/e^β[ε_v( k_2) - μ_v]+1 - 1/e^β[ε_c( k_1) - μ_c]+1) =  n_v( k_2)-n_c( k_1)/ε_c( k_1) - ε_v( k_2) -μ - i Ω ,

as given in Eqs. (<ref>)–(<ref>) of the main text.



Analogously we perform the calculation of the loop for the fourth-order term, presented in Fig. <ref>(b):

    ℒ_ k_1 k_2 k_3 k_4^Ω_1Ω_2Ω_3 =T ∑_ω'𝒢_0( k_1, ω^')_11𝒢_0( k_2, ω^'- Ω_1)_22𝒢_0( k_3, ω^'- Ω_1 + Ω_2)_11𝒢_0( k_4, ω^' - Ω_1 + Ω_2 - Ω_3)_22
    
        = T1/2π i1/T∮ dz 1/e^β z+11/[-z  + ε_c( k_1)- μ_c][-z + i Ω_1 + ε_v( k_2)- μ_v]
    ×1/[-z + i Ω_1 - iΩ_2 + ε_c( k_3)- μ_c][-z + i Ω_1 - iΩ_2 + i Ω_3 +ε_v( k_4)- μ_v]
    
        =-{n_c( k_1)/[iΩ_1 + ε_v( k_2)-ε_c( k_1)+ μ][iΩ_1- i Ω_2 + ε_c( k_3)-ε_c( k_1)][iΩ_1- i Ω_2 +i Ω_3 + ε_v( k_4)-ε_c( k_1) +μ].
    
        +n_v( k_2)/[-iΩ_1 + ε_c( k_1)-ε_v( k_2) - μ][- i Ω_2 + ε_c( k_3)-ε_v( k_2)-μ][- i Ω_2 +i Ω_3 + ε_v( k_4)-ε_v( k_2)]
    
        + n_c( k_3)/[-iΩ_1 + iΩ_2+ ε_c( k_1)-ε_c( k_3)][i Ω_2 + ε_v( k_2)-ε_c( k_3)+ μ][i Ω_3 + ε_v( k_4)-ε_c( k_3)+μ]
    .+n_v( k_4)/[-iΩ_1 + iΩ_2- iΩ_3+ ε_c( k_1)-ε_v( k_4)-μ][i Ω_2 -iΩ_3+ ε_v( k_2)-ε_v( k_4)][-i Ω_3 + ε_c( k_3)-ε_v( k_4)-μ]}
    
        =-1/iΩ_1- i Ω_2 + ε_c( k_3)-ε_c( k_1){n_c( k_1)/[iΩ_1 + ε_v( k_2)-ε_c( k_1)+μ][iΩ_1- i Ω_2 +i Ω_3 + ε_v( k_4)-ε_c( k_1)+μ].
        
    .-n_c( k_3)/[i Ω_2 + ε_v( k_2)-ε_c( k_3)+μ][i Ω_3 + ε_v( k_4)-ε_c( k_3)+μ]}
    
    
    
        - 1/- i Ω_2 +i Ω_3 + ε_v( k_4)-ε_v( k_2){n_v( k_2)/[-iΩ_1 + ε_c( k_1)-ε_v( k_2)-μ][- i Ω_2 + ε_c( k_3)-ε_v( k_2)-μ]. 
    
        -. n_v( k_4)/[-iΩ_1 + iΩ_2- iΩ_3+ ε_c( k_1)-ε_v( k_4)-μ][-i Ω_3 + ε_c( k_3)-ε_v( k_4)-μ]}.

Introducing for convenience the notation

    B_k_1 k_2   = -iΩ_1 + ε_c( k_1) - ε_v( k_2)-μ,     
        B_k_3 k_2 = -iΩ_2 + ε_c( k_3) - ε_v( k_2) -μ,
    
        B_k_3 k_4   = -iΩ_3 + ε_c( k_3) - ε_v( k_4) -μ,     B_k_1 k_4 = -iΩ_1 + -iΩ_2 -iΩ_3 + ε_c( k_1) - ε_v( k_4)-μ,

we express the fourth-order loop in the following form:

    ℒ_ k_1 k_2 k_3 k_4^Ω_1Ω_2Ω_3 = 1/B_k_1k_2- B_k_3 k_2[n_c( k_1)/B_k_1k_2 B_k_1k_4 - n_c( k_3)/B_k_3k_2B_k_3k_4] -1/B_k_3k_2 - B_k_3k_4[n_v( k_2)/B_k_1k_2B_k_3k_2 - n_v( k_4)/B_k_1k_4B_k_3k_4].



The third-order loops that are of interest to our calculations presented in Appendix <ref> are shown in Fig. <ref>(c) and are defined as

    [Π_121]_ k_1 k_2 k_3^ Ω_1  Ω_2   =    T∑_ω^'𝒢_0( k_1, ω^')_11𝒢_0( k_2,ω^'-Ω_1)_22𝒢_0( k_3, ω^'- Ω_1 + Ω_2)_11,  
    [Π_122]_ k_1 k_2 k_3^ Ω_1  Ω_2   =    T∑_ω^'𝒢_0( k_1, ω^')_11𝒢_0( k_2, ω^'- Ω_1)_22𝒢_0( k_3, ω^'- Ω_2)_22.

The calculation is performed in the same way:

    [Π_121]_ k_1 k_2 k_3^ Ω_1  Ω_2 = T1/2 π i1/T∮  dz 1/e^β z+1 1/-z +ε_c( k_1) -μ_c 1/-z + i Ω_1 +ε_v( k_2) -μ_v 1/-z + i Ω_1 - i Ω_2  +ε_c( k_3) -μ_c
    
         = - {n_c( k_1)/[i Ω_1 + ε_v( k_2) - ε_c( k_1) + μ][i Ω_1 - i Ω_2 +ε_c( k_3) - ε_c( k_1)] +  . 
    . n_v( k_2)/[-iΩ_1 +ε_c( k_1) -ε_v( k_2) -μ][-iΩ_2 +ε_c( k_3) -ε_v( k_2) -μ] + n_c( k_3)/[-i Ω_1 + i Ω_2 +ε_c( k_1) -ε_c( k_3)][i Ω_2 +ε_v( k_2) -ε_c( k_3) +μ]}


    = n_c( k_1)/B_k_1k_2(B_k_1k_2 - B_k_3k_2) + n_v( k_2)/B_k_1k_2 B_k_3k_2 + n_c( k_3)/B_k_3k_2(B_k_3k_2-B_k_1k_2),


    [Π_122]_ k_1 k_2 k_3^ Ω_1  Ω_2 = - {n_c( k_1)/[i Ω_1 + ε_v( k_2) -ε_c( k_1) + μ][i Ω_1 - iΩ_2 + ε_v( k_3) - ε_c( k_1)+μ] +.
    .  n_v( k_2)/[-i Ω_1 +ε_c( k_1) -ε_v( k_2) -μ][-i Ω_2 +ε_v( k_3) -ε_v( k_2)] + n_v( k_3)/[-i Ω_1 + iΩ_2 +ε_c( k_1) -ε_v( k_3) -μ][i Ω_2 +ε_v( k_2) -ε_v( k_3)]}
    
        = n_c( k_1)/B_k_1k_2B_k_1k_3 +n_v( k_2)/B_k_1k_2(B_k_1k_2-B_k_1k_3) +n_v( k_3)/B_k_1k_3(B_k_1k_3-B_k_1k_2).




§ CALCULATIONS AT ZERO TEMPERATURE

This Appendix is devoted to simplification of the general expressions obtained in Sec. <ref> for the case T=0 which allows to obtain some analytical results comparable with those exciting in literature (see discussion in the main text). 

Assuming the limit of a large system size, the Fourier transforms can be rewritten as

    
    Ψ   _i( r, τ) =∫d k/(2π)^2dω/2π Ψ_i( k, ω) e^i  k· x - iωτ,
       Δ   (x, y) = ∫d k/(2π)^2d k^'/(2π)^2dΩ/2π Δ( k,  k^', Ω) e^i k· x - i  k^'· y - iΩτ,
       ϕ   (x, y) =∫d k/(2π)^2d k^'/(2π)^2dΩ/2π ϕ( k,  k^', Ω) e^i k· x - i  k^'· y - iΩτ.


At zero temperature n_v=1, n_c=0, and one can derive from the expansion series the analytical expressions for all the loops considered in Appendix <ref>,

    
    𝒜_ k_1 k_2^  Ω   =   1/ε_c( k_1) - ε_v( k_2) - μ- i Ω, 
    [Π_121]_ k_1 k_2 k_3^ Ω_1  Ω_2   =   1/[-i Ω_1 + ε_c( k_1) - ε_v( k_2) - μ][-i Ω_2 + ε_c( k_3) - ε_v( k_2) - μ], 
    [Π_122]_ k_1 k_2 k_3^ Ω_1  Ω_2   =   1/[-i Ω_1 + ε_c( k_1) - ε_v( k_2) - μ][-i Ω_1 + i Ω_2 + ε_c( k_1) - ε_v( k_3) - μ].
    ℒ_ k_1 k_2 k_3 k_4^Ω_1Ω_2Ω_3   =   -i Ω_1 +i Ω_2 - i Ω_3 + ε_c( k_1) - ε_v( k_4) -i Ω_3 + ε_c( k_3) - ε_v( k_4)  - 2μ/[-i Ω_1 + ε_c( k_1) - ε_v( k_2) - μ][-i Ω_2 + ε_c( k_3) - ε_v( k_2) - μ]
              ×1/[-i Ω_3 + ε_c( k_3) - ε_v( k_4) - μ][-i Ω_1 +i Ω_2 + i Ω_3 + ε_c( k_4) - ε_v( k_4) - μ].


While the fourth–order loop formula (<ref>) is still quite sophisticated, within the saddle-point approximation all the denominators vanish and we get:

    ϕ( k_1, k_2,Ω) = -i/𝒜_ k_1 k_2^  Ω Δ( k_1, k_2,Ω).


Going back to Eq. (<ref>), we turn to the position basis

for the relative motion r and write down the generalized Wannier equation in terms of the exciton field Δ( r,  k, Ω):

    [-i Ω + E_ g + ħ^2  k^2/2 M - ħ^2/2m ∇^2 - V(r)  - μ]Δ( r,  k, Ω) -g_ R^2Δ( r=0,  k, Ω)/E_ ph( k) -μ -i Ω=0,

where m the the reduced mass and M is the exciton mass.
It is worth noting that the same approach has been applied in Ref. <cit.>, albeit for the one-mode fields within the semiconductor Bloch equations approach. The field Δ( r,  k, Ω) in such mixed representation can be 
regarded as g_ R^2Δ( r=0,  k, Ω)/[E_ ph( k) -μ -i Ω] G_ k, Ω( r,  r^'), where the Green's function is defined from the equation:

    [-i Ω + E_ g  - μ+ ħ^2  k^2/2 M - ħ^2/2m∇^2 - V(r)]G_ k, Ω( r,  r^') = δ( r-r^').

Here we suggest the solution in terms of the ν s–state of the hydrogen-like equation, although 
the potential of interaction for χ_ν( r) can be considered to be of an arbitrary shape: [- ħ^2∇^2/2m - V(r)]χ_ν( r) = E_νχ_ν( r). The Green's function then has the form:

    G_ k, Ω( r,  r^') = ∑_ν A^(ν)_ k, Ω( r')χ_ν( r) = ∑_νχ^*_ν( r^') χ_ν( r)/-iΩ + ħ^2 k^2/2M + E_ g - μ + E_ν,

and the exciton field turns into

    Δ( r,  k, Ω) = g_ R^2Δ( r=0,  k, Ω)/E_ ph( k) -μ -i Ω  G_ k, Ω( r, 0).

As one 
sees, the sum in the derived expression (<ref>) diverges in the point r=0 as ∼ -lnr, as becomes clear from the form of Eq. (<ref>) when estimating the trivial limit of k =0, Ω=0 at r→0, when one can neglect the potential compared to the delta-function in the right-hand side.
This divergence is known <cit.> and arises because 
the attractive potential of the delta-functional form produced by the photon field leads to the reduction of the mean electron-hole separation (see e.g. <cit.>).

The exciton-photon action in momentum 
basis for the c.m. motion 
takes the form:

    𝒮= ∑_ν∫d  k/(2π)^2∫d Ω/2 π{Δ̅( r=0,  k, Ω)
        |g_ R^2/E_ ph( k)- μ -i Ω|^2 |χ_ν( r=0)|^2/iΩ + E_ g - μ + ħ^2  k^2/2M + E_ν Δ( r=0,  k, Ω)   .
    . + g_ R[ Ψ_ ph( k, Ω) g_ R^2 |χ_ν( r=0)|^2/E_ ph( k) -μ- i Ω Δ( r = 0,  k, Ω)/-iΩ+ E_ g-μ+ħ^2  k^2/2M + E_ν+Δ̅( r =0,  k, Ω)/iΩ+ E_ g-μ+ħ^2 k^2/2M + E_ν g_ R^2  |χ_ν( r = 0)|^2/E_ ph( k) -μ+ i ΩΨ_ ph( k, Ω) ]}
     + ∫d  k/(2π)^2∫d Ω/2 π  Ψ_ ph( k,Ω)[- i Ω + E_ ph( k) -μ] Ψ_ ph( k,Ω)
    
        + 1/2∫d p/(2π)^2∫d l_1 d l_2 d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3[ -i Ω_1 + i Ω_2 - i Ω_3 + ε_c( p+ l_1/2) - ε_v( p+ l_2 - l_3 - l_1/2) - i Ω_3 + ε_c( p+ l_2- l_1/2) .
    .- ε_v( p+ l_2- l_3 - l_1/2)]  Δ( p,  l_1, Ω_1)Δ̅(  p+ l_2 - l_1/2,  l_2, Ω_2)
        Δ( p+ l_2 - l_3+ l_1/2,  l_3, Ω_3) Δ̅(  p+ l_2 - l_3/2,  l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3).

As one can see from three first terms in Eq. (<ref>), it is natural to introduce 
a field corresponding to the exciton as a whole (i.e. depending only on the total momentum):

    C^(ν)( k,Ω) =  g_ R^2 χ^*_ν( r =0)/E_ ph( k)- μ -i ΩΔ( r=0,  k, Ω)/[-iΩ + E_ g  - μ + ħ^2 k^2/2M + E_ν],

hence the exciton field Δ can be rewritten as

    Δ( p,  k, Ω) = ∑_ν C^(ν)( k, Ω) χ_ν( p).

Provided all these simplifications and with the exciton interaction term taken in the saddle point according to Eq. (<ref>), the action (<ref>) takes the form 

    𝒮= ∑_ν∫d  k/(2π)^2∫d Ω/2 π{C^(ν)( k, Ω)(-iΩ+ E_ g-μ+ħ^2 k^2/2 M+ E_ν)C^(ν)( k, Ω) + g_ R[C^(ν)( k, Ω)χ^*_ν( r=0)  Ψ_ ph( k, Ω) ..
         
    .. + Ψ_ ph( k, Ω) C^(ν)( k, Ω) χ_ν( r=0) ]} + ∫d  k/(2π)^2∫d Ω/2 π  Ψ_ ph( k,Ω)[- i Ω + E_ ph( k) -μ]Ψ_ ph( k,Ω)
         
    
         +1/2∑_ν_1…ν_4∫d p/(2π)^2∫d l_1 d l_2 d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3 C^(ν_1)( l_1, Ω_1) χ_ν_1( p)C^(ν_2)( l_2, Ω_2) χ^*_ν_2( p +  l_2 - l_1/2)C^(ν_3)( l_3, Ω_3) χ_ν_3( p +  l_2 -  l_3 + l_1/2) 
    ×
        2 ∫d q/(2π)^2  V( q)  C^(ν_4)( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3 )χ^*_ν_4( p+ l_2 - l_3/2- q)


    -1/2∑_ν_1…ν_3∫d p/(2π)^2∫d l_1 d l_2 d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3 [ C^(ν_1)( l_1, Ω_1) χ_ν_1( p)C^(ν_2)( l_2, Ω_2) χ^*_ν_2( p+ l_2 - l_1/2) . 
    . ×  C^(ν_3)( l_3, Ω_3) χ_ν_3( p+ l_2 - l_3 + l_1/2) g_ RΨ_ ph( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3 ) + c.c.],

where the last term describes the so-called “saturation” or, 
to be precise, the exciton-assisted exciton-photon coupling. 
We note here that to the final expression for action one should add the contribution to the exciton-exciton interaction arising from the electron screening derived in Appendix <ref>:

    𝒮^'_ scr = - 1/2∑_ν_1…ν_4∫d p/(2π)^2∫d l_1 d l_2 d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3  C^(ν_1)( l_1, Ω_1)χ_ν_1( p)C^(ν_2)( l_2, Ω_2) χ^*_ν_2( p+ l_2 - l_1/2) 
    × 2 ∫d q/(2π)^2  V( q)C^(ν_3)( l_3, Ω_3) χ_ν_3( p - q+ l_2 - l_3 + l_1/2) C^(ν_4)( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3)χ^*_ν_4( p - q+ l_2- l_3/2).

The action obtained as a sum of Eqs. (<ref>) and (<ref>) has the form quite similar to the commonly-used action, but containing the sums over the exciton s–states.




§ THE “RIGID” EXCITON LIMIT


In the following we restrict ourselves with the standard case in which only the 1s–exciton state is assumed to appear in the system, revealing how the expressions for interactions can be transformed into the widely known values of interaction constants.  
To this end, it is convenient to restrict ourselves in Eq. (<ref>) to χ^(1)_ k,Ω( p)≡χ_ k( p) and the strong-coupling regime ħΩ_ R≪ E_ b, and thus neglect the term ∼(ħΩ_ R)^2/E_ ph. In this case, without the van-der-Waals term (see Appendix <ref>) Eq. (<ref>) turns into:

    (E_g +ħ^2  k^2/2 M +ħ^2  p^2 /2m - μ - i Ω)Δ( p, k, Ω) -[n_v( p - m/m_v k) - n_c( p + m/m_c k)]∫d q/(2π)^2V( p - q) Δ(  q, k,Ω) =0

and can be solved using the ansatz Δ(  p, k,Ω) = χ_ k( p)C( k, Ω) where the 1s exciton wavefunction χ_ k( p) is modified by temperature.

At T=0, the saddle-point equation for configuration of the field Δ turns into the standard Wannier equation:

    (E_ g +ħ^2  k^2/2 M +ħ^2  p^2 /2m - μ - i Ω)Δ( p, k, Ω) -∫d q/(2π)^2V( p - q) Δ(  q, k,Ω) =0.

The solution is  Δ( p,  k, Ω) = χ( p)C( k, Ω) where χ( p) is the 1s–exciton  wavefunction and the excited s–states are neglected:

    ( p^2/2m - E_ 1s)χ( p) = ∫d q/(2π)^2V( p - q)χ( q).

In this case, the saddle-point equations (<ref>) and (<ref>) reduce to the standard Hopfield equations, and  the exciton-photon action takes the well-known form <cit.>:

    𝒮= ∫d  k/(2π)^2∫d Ω/2 π{C( k, Ω) (-iΩ + E_ g  - μ + ħ^2  k^2/2 M + E_ 1s)C( k, Ω)+  Ψ_ ph( k,Ω)(- i Ω + E_ ph( k) -μ)Ψ_ ph( k,Ω) .
    .+ħΩ_ R/2[C( k, Ω)Ψ_ ph( k, Ω)  + Ψ_ ph( k, Ω) C( k, Ω) ] }
    
         + 1/2∫d l_1 d l_2 d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3{V_ ex^tot( l_1, l_2,  l_3)  C( l_1, Ω_1) C( l_2, Ω_2) C( l_3, Ω_3)C( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3) .
    -
         . [V_ sat( l_1,  l_2,  l_3)C( l_1, Ω_1) C( l_2, Ω_2) C( l_3, Ω_3)Ψ_ ph( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3)  + c.c.]}

with ħΩ_ R/2 = g_ R∑_ qχ( q) = g_ Rχ( r=0) being the experimentally relevant Rabi-splitting (whereas g_ R is the bare electron-hole-photon coupling rate), and the notations

    V^tot_ ex( l_1, l_2, l_3)    ≡ 2 ∫d p/(2π)^2∫d q/(2π)^2 V( q) χ( p) χ̅( p+ l_2 - l_1/2)χ̅( p+ l_2 - l_3/2- q)[χ( p+ l_2 - l_3 + l_1/2) -χ( p+ l_2 - l_3 + l_1/2- q)]
         
    
        V_ sat( l_1, l_2, l_3)    ≡ g_ R∫d p/(2π)^2 χ( p)χ̅( p+ l_2 - l_1/2) χ( p+ l_2 - l_3 + l_1/2).

One sees that depending on the material and the interaction potential between the charges within the exciton, the shape of χ( p) may vary thus alternating both the exciton-exciton (interaction) and the saturation nonlinearities in the polariton system. 


We note as well that the exciton-exciton interaction term including Eq. (<ref>) is usually calculated 
in the operator formalism <cit.>. Here we stress that within our approach we derive excitonic rather than polaritonic interactions (if one instead solves the Bethe-Salpeter equation modified by light (<ref>), the obtained interactions of the renormalized field would be “polaritonic”).




Finally, assuming that 
the c.m. momenta l_i (i=1,2,3) are negligibly small compared to the momentum of internal exciton motion, the  exciton-exciton interaction constant obtained from Eq. (<ref>) coincides with the well-known result:

    g_ ex= 2 ∫d p/(2π)^2∫d q/(2π)^2 V( q)χ( p)χ^*( p)χ^*( p - q)[χ( p )-χ( p -q )].






The 
saturation interaction constant is obtained similarly from Eq. (<ref>): 

    g_ sat = ħΩ_ R/2χ( r=0)∫d p/(2π)^2 χ( p)χ^*( p)χ( p),

which coincides with the expressions derived using alternative approaches in Refs. <cit.>.


Starting from the action (<ref>) 
and considering 
the limit of Eqs. (<ref>), (<ref>),

    𝒮 = ∫d  k/(2π)^2∫d Ω/2 π{C( k, Ω)(-iΩ + E_ g  - μ + ħ^2  k^2/2 M + E_ 1s)C( k, Ω) +  Ψ_ ph( k,Ω)[- i Ω + E_ ph( k) -μ]Ψ_ ph( k,Ω) +.
    . + ħΩ_ R/2[C( k, Ω)Ψ_ ph( k, Ω)  + Ψ_ ph( k, Ω) C( k, Ω) ] }
         + ∫d l_1d l_2d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3{g_ ex/2C( l_1, Ω_1) C( l_2, Ω_2)  C( l_3, Ω_3) . 
    ×C( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3) . -  g_ sat/2
         [ C( l_1, Ω_1) C( l_2, Ω_2)  C( l_3, Ω_3)Ψ_ ph( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3 ) +  c.c.]},

we perform the saddle-point approximation assuming the system to be in equilibrium iΩ = 0:

    
    . δ𝒮_ eff/δC|_C=Ψ_ ex
    Ψ_ ph     =    (E_ g + E_1s-μ) Ψ_ ex + ħΩ_ R/2Ψ_ ph +g_ ex |Ψ_ ex|^2 Ψ_ ex -g_ sat |Ψ_ ex|^2 Ψ_ ph - g_ sat/2(Ψ_ ex)^2 Ψ_ ph= 0, 
    . δ𝒮_ eff/Ψ_ ph^0|_C=Ψ_ ex
    Ψ_ ph     =    (E_ ph-μ) Ψ_ ph +ħΩ_ R/2Ψ_ ex - g_ sat/2 |Ψ_ ex|^2 Ψ_ ex=0


where Ψ_ ex(ph) are the equilibrium exciton (photon) fields. 
The relation between the two 
can be defined as:

    Ψ_ ph = Ψ_ ex/μ  - E_ ph^0(ħΩ_ R/2 - g_ sat|Ψ_ ex|^2/2),

and the phase difference between 
them being 0 or π depending on μ. Solving these algebraic equations, one can calculate the bottoms of dispersion laws of the two new normal modes, i.e. lower (LP) and upper (UP) polaritons, which depend on the uniform exciton density n_ ex^0:

    μ_ LP(UP)(n_ ex^0) = E_ g + E_1s + E_ ph^0+ g_ ex n_ ex^0/2∓1/2√(( g_ ex n_ ex^0 - Δ)^2 + (ħΩ_ R - g_ sat n_ ex^0)(ħΩ_ R - 3 g_ sat n_ ex^0)),

where E_ ph^0 - (E_ g + E_1s)= Δ is usually regarded as the (constant) energy detuning between photons and excitons at k=0.

Rewriting the formulae in terms of the renormalized (density-dependent) detuning Δ→Δ(n_ ex^0) = Δ - g_ ex n_ ex^0 and Rabi-splitting ħΩ_R →ħΩ_R(n_ ex^0) = √((ħΩ_ R - g_ sat n_ ex^0)(ħΩ_ R - 3 g_ sat n_ ex^0)), we get:

    μ_ LP(UP)(n_ ex^0) = E_ g + E_1s+ g_ ex n_ ex^0 + Δ (n_ ex^0)/2∓1/2√(Δ^2(n_ ex^0) + [ħΩ_R(n_ ex^0)]^2).

One sees that the separation of the bottoms of dispersions (at k=0) shifts with the exciton density as

    μ_ UP(n_ ex^0)- μ_ LP(n_ ex^0) = √(Δ^2(n_ ex^0) + [ħΩ_R(n_ ex^0)]^2).

At the same time, from the experimental point of view a more relevant quantity characterizing the exciton-photon conversion rate is the distance between the polariton branches at the anticrossing point (which is shifting itself with the change of the exciton density, since the detuning is also changing). It can be defined as follows: 
    min[E_ UP( k)- E_ LP( k)] = √(Δ_ k^2(n_ ex^0) + [ħΩ_R(n_ ex^0)]^2)  with Δ_ k(n_ ex^0) = Δ - g_ ex n_ ex^0 + ħ^2  k^2/2 m_ ph - ħ^2  k^2/2 m_ ex.




§ CONSIDERATION OF SPINS

For purely excitonic systems, the treatment of spins within functional integration approach was briefly mentioned in Ref. <cit.> without any detailed discussion. The aim of this Appendix is to study the influence of spins on the effective exciton-photon action in polariton systems.
We restrict ourselves to the simplest case, considering T=0 and rigid  1s–excitons. Moreover, we will focus on the corrections to the expressions derived above arising due to dark exciton states, so the density channels of pairing addressed in Appendix <ref> at this stage are neglected for convenience. Below we shortly discuss the screening of exciton interaction for excitons with different spins. 
Even such a simple case upon examination turns out to be intricate.

Our starting point is the action without spin-orbit coupling:

    𝒮[Ψ_c, Ψ_v, Ψ_ ph] = ∫d r∫_0^β dτ{[ Ψ_c ↑(x);  Ψ_c↓(x);  Ψ_v↑(x); Ψ_v ↓(x) ]^ T[ ∂_τ + ε_c ↑(k̂)               0               0               0;               0  ∂_τ + ε_c↓(k̂)               0               0;               0               0  ∂_τ + ε_v↑(k̂)               0;               0               0               0 ∂_τ + ε_v ↓(k̂) ][ Ψ_c ↑(x); Ψ_c ↓(x); Ψ_v ↑(x); Ψ_v ↓(x) ].
            
    
            . + g_ R[Ψ_ ph ↑(x)Ψ_c ↓(x)Ψ_v ↑(x) +Ψ_c ↓(x)Ψ_v ↑(x)Ψ_ ph↑(x) + Ψ_ph ↓(x)Ψ_c ↑(x)Ψ_v ↓(x) +Ψ_c ↑(x)Ψ_v ↓(x)Ψ_ ph↓(x) ]}
    
            +  ∑_σ∫ d r∫_0^β dτ Ψ_ ph, σ(x)[∂_τ + E_ ph σ(k̂)] Ψ_ ph, σ(x)
            +1/2∑_i, j∑_σ, σ^'∫ d rd r^'∫_0^β dτ dτ^' V(x-x^') Ψ_i σ(x)Ψ_i σ(x) Ψ_j σ^'(x^')Ψ_j σ^'(x^'),

where σ = ↑, ↓ is a spin index, x=( r,t) and i,j=c,v as before, ε_c(v) ↑ (↓)( k) = ± E_ g/2 ±ħ^2  k^2/2m_c(v)↑ (↓) - μ_c(v)↑ (↓), E_ ph ↑(↓)( k) = E_ ph^0 + ħ^2  k^2/2 m_ ph -μ_ ph ↑ (↓). Here, in general case, we assume the chemical potentials and effective masses of electrons with different spins in conduction and valence bands being different.
We introduce four exciton fields corresponding to the bright and dark excitons. For TMD materials, where the electrons in both valence and conduction bands have spin projections ±1/2, one has:
    
    Δ_+1(x, y) = Ψ_v↑ (y)Ψ_c ↑(x), Δ_-1(x, y) = Ψ_v↓ (y)Ψ_c ↓(x), Δ_0(x, y) = Ψ_v↑ (y)Ψ_c ↓(x),  Δ_0̅(x, y) = Ψ_v↓ (y)Ψ_c ↑(x),

while for GaAs, where we consider only heavy holes, i.e. electrons in the valence band with spin projections ± 3/2 which 
correspond to Ψ_v ↑ and Ψ_v ↓, respectively,
    
    Δ_+1(x, y) = Ψ_v↑ (y)Ψ_c ↓(x), Δ_-1(x, y) = Ψ_v↓ (y)Ψ_c ↑(x), Δ_+2(x, y) = Ψ_v↑ (y)Ψ_c ↑(x),  Δ_-2(x, y) = Ψ_v↓ (y)Ψ_c ↓(x).

In the following, we will not make a distinction between the materials and will generically denote the two dark fields as Δ_ d,d̅. Along with the four exciton fields and their conjugates, in the same way as in the spinless case (using 8 δ–functions, see the main text), we have 4 auxiliary fields ϕ_±1, ϕ_ d,d̅ and their conjugates. 
The action (<ref>) takes the form:

    𝒮[Ψ_c,v, Ψ_ ph, ϕ_±1, d,d̅, Δ_±1, d,d̅] = ∑_ k,ω[  Ψ_c ↑(k);  Ψ_c ↓(k);  Ψ_v ↑(k); Ψ_v  ↓(k) ]^T[  -i ω + ε_c ↑(k)                0                0                0;                0    -i ω+ ε_c↓(k)                0                0;                0                0  -i ω + ε_v↑( k)                0;                0                0                0 -i ω + ε_v ↓( k) ][ Ψ_c  ↑(k); Ψ_c  ↓(k); Ψ_v  ↑(k); Ψ_v  ↓(k) ]
    
        + i√(T)∑_ k_1, k_2 
    ω_1, ω_2[ Ψ_c ↑(k_1); Ψ_c ↓(k_1); Ψ_v ↑(k_1); Ψ_v ↓(k_1) ]^Tδ𝒢^-1(k_1, k_2) [ Ψ_c  ↑(k_2); Ψ_c  ↓(k_2); Ψ_v  ↑(k_2); Ψ_v  ↓(k_2) ]
            -i ∑_σ = ±1, 
     d,d̅∑_ k_1, k_2 
    Ω[ ϕ̅_σ( k_1, k_2, Ω)Δ_σ( k_1, k_2, Ω) + Δ̅_σ( k_1, k_2, Ω) ϕ_σ( k_1, k_2, Ω) ] 
            
    
            - ∑_σ = ±1, 
     d,d̅∑_ k_1… k_4 
    Ω V( k_1 -  k_2) Δ̅_σ( k_1, k_4,Ω) Δ_σ( k_2, k_3,Ω) δ( k_1+ k_3,  k_2+ k_4) + ∑_σ = ±1∑_ k, ΩΨ_ ph,σ(k)[- i Ω + E_ ph( k) -μ]Ψ_ ph,σ(k)  
            
    
            + g_ R/√(S)∑_σ = ±1∑_ k_1, k_2 
    Ω[Ψ_ ph,σ( k_1- k_2,Ω) Δ_σ( k_1, k_2,Ω) + Δ̅_σ( k_1, k_2,Ω)Ψ_ ph, σ( k_1- k_2,Ω)],

where

    δ𝒢^-1(k_1, k_2)= [                           0                           0   ϕ_ d ( k_1, k_2, ω_1-ω_2)  ϕ_-1 ( k_1, k_2,  ω_1-ω_2);                           0                           0  ϕ_+1 ( k_1, k_2,  ω_1-ω_2)  ϕ_d̅ ( k_1, k_2,  ω_1-ω_2); ϕ̅_ d ( k_2, k_1,  ω_2-ω_1)  ϕ̅_+1 ( k_2, k_1, ω_2-ω_1)                           0                           0;  ϕ̅_-1 ( k_2, k_1, ω_2-ω_1)  ϕ̅_d̅ ( k_2, k_1, ω_2-ω_1)                           0                           0 ].

In the above expressions, we used the spin composition for GaAs, but the result is similarly applicable to TMD materials and can be obtained by replacing ϕ_ d, d̅→ϕ_± 1 and ϕ_± 1→ϕ_0, 0̅ in (<ref>).
After integrating over the fermionic fields, expanding Tr ln𝒢^-1 up to the second order in δ𝒢^-1 and performing the saddle-point approximation, we obtain the relation between ϕ_σ and Δ_σ which is the same as for the previously considered case, but containing spin indices:

    ϕ_σ( k_1, k_2,Ω) = Δ_σ( k_1, k_2,Ω)/i[𝒜_σ]_ k_1 k_2^  Ω,   ϕ̅_σ( k_1, k_2,Ω) = Δ̅_σ( k_1, k_2,Ω)/i[𝒜_σ]_ k_1 k_2^  Ω,

where at T=0

    1/[𝒜_σ]_ k_1 k_2^  Ω = 
            -i Ω+ ε_c↑( k_1) - ε_v↑( k_2)   if σ =  d
    
             -i Ω +ε_c↓( k_1) - ε_v↑( k_2)   if σ = +1 
    
             -i Ω +ε_c↑( k_1) - ε_v↓( k_2)   if σ = -1 
    
              -i Ω +ε_c↓( k_1) - ε_v↓( k_2)   if σ = d̅.

In principle, the expressions in the right-hand side of Eq. (<ref>) are the exciton dispersions. The chemical potentials of the exciton fields with different spin projections read: μ_c↑ - μ_v ↑ = μ_ d, μ_c↓ - μ_v ↑=μ_+1, μ_c↑ - μ_v ↓ =μ_-1, μ_c↓ - μ_v ↓ = μ_d̅. Moreover, under the assumption of the thermodynamic equilibrium between photons and excitons, one gets μ_+1 = μ_ ph ↑, μ_-1 =μ_ ph ↓. Then the saddle-point approximation for exciton fields yields:


    Δ_± 1( p, k, Ω)/[𝒜_±1]_ k_1 k_2^  Ω    + g_ R/√(S)Ψ_ ph, ± 1( k, Ω) -∑_ qV( p - q) Δ_± 1(  q, k,Ω)=0, 
       Δ_ d(d̅)( p, k, Ω)/[𝒜_ d(d̅)]_ k_1 k_2^  Ω    - ∑_ qV( p - q) Δ_ d(d̅)(  q, k,Ω)=0,

 where p and k are the relative motion of electrons and the exciton c.m. momenta, respectively. In the c.m. frame the second-order loop 𝒜^-1 takes the form:

    1/[𝒜_σ]_ k_1 k_2^  Ω = E_ g -μ_σ + ħ^2  p^2/2 m_σ + ħ^2  k^2/2 M_σ,

with m_σ being the reduced mass of the σ-exciton field and M_σ the σ-exciton mass. 

In the rigid exciton limit the variables can be divided Δ_σ( p,  k, Ω) = χ_σ( p)C_σ( k, Ω), similarly to the spinless case. The 1s–state wavefunctions χ_σ( p) for different fields are determined by the band splitting: 

    p^2/2m_σχ_σ( p) - ∑_ qV( p-q)χ_σ( q)= E_ 1s^σχ_σ( p).

If the bands for different spins coincide, the expressions are trivial since 𝒜 becomes independent of spin.

In principle, the second-order expansion does not deviate significantly from the spinless case and all the loop calculations are straightforward and similar to those presented in Appendix <ref>, albeit cumbersome. Here we present only the final result. 

For the fourth-order term in the expansion series (1/4) Tr(𝒢_0 δ𝒢^-1)^4
one has 32 loops. Four 
of them describe bright-exciton interactions (two for |ϕ_+1|^4 and two for |ϕ_-1|^4):

    1/4 Tr(𝒢_0δ𝒢^-1)^4 =2/4(i√(T))^4 ∑_ k_1… k_4∑_Ω_1…Ω_3ϕ_± 1( k_1, k_2, Ω_1)ϕ̅_± 1( k_3,  k_2, Ω_2) ϕ_± 1( k_3,  k_4, Ω_3)ϕ̅_± 1( k_1,  k_4, Ω_1 + Ω_3 -Ω_2) 
    ×∑_ω^'𝒢_0( k_1, ω^')_22(11)𝒢_0( k_2, ω^' - Ω_1)_33(44)𝒢_0( k_3, ω^'- Ω_1+ Ω_2)_22(11)𝒢_0( k_4, ω^'- Ω_1+ Ω_2- Ω_3)_33(44)
    
       = 1/4T ∑_ k_1… k_4∑_Ω_1…Ω_3Δ_± 1( k_1, k_2, Ω_1)Δ̅_± 1( k_3,  k_2, Ω_2) Δ_± 1( k_3,  k_4, Ω_3)Δ̅_± 1( k_1,  k_4, Ω_1 + Ω_3 -Ω_2) ×
    [-i Ω_1 + ε_c↓(↑)( k_1) - ε_v↑(↓)( k_2) - i Ω_3 + ε_c↓(↑)( k_3) - ε_v↑(↓)( k_4) - i Ω_2 + ε_c↓(↑)( k_3) - ε_v↑(↓)( k_2) 
       - i Ω_1 + iΩ_2 - i Ω_3 + ε_c↓(↑)( k_1) - ε_v↑(↓)( k_4)]

which should be treated according to the saddle-point equation leading to the exciton-exciton interaction and saturation terms, similarly with the formulae from Appendix <ref> (it is worth noting that the obtained expression has symmetric form, for convenience).

Dark-exciton interactions |ϕ_ d,d̅|^4 contribution has generally the same symmetrized shape (4 diagrams):

    2/4T ∑_ k_1… k_4∑_Ω_1…Ω_3Δ_ d,d̅( k_1, k_2, Ω_1)Δ̅_ d,d̅( k_3,  k_2, Ω_2) Δ_ d,d̅( k_3,  k_4, Ω_3)Δ̅_ d,d̅( k_1,  k_4, Ω_1 + Ω_3 -Ω_2) 
    ×[-i Ω_1 + ε_c↑(↓)( k_1) - ε_v↑(↓)( k_2)  -i Ω_2 + ε_c↑(↓)( k_3) - ε_v↑(↓)( k_2)].

Due to the form of the saddle-point equation, this contribution is not leading to saturation,  which is anticipated since the dark excitons are not coupled to photons.

Interestingly, bright-bright and dark-dark interactions do not involve excitons with opposite spins, as long as the Hamiltonian does not 
contain the spin-orbit coupling term [i.e. the Green's function 𝒢^-1_0 needs to contain non-zero (𝒢^-1_0)_13 or(and)  (𝒢^-1_0)_24 elements]. In existing literature <cit.> such interactions are also attributed to the biexciton resonance (excitons with opposite spins forming a biecxiton), which cannot be treated within our approach since it does not consider multiple-exciton bound states.

The goal of this Appendix is however to take into account the dark states when considering bright-dark exciton interactions and especially saturation, hence the main subject of our consideration are the terms |ϕ_+1|^2 |ϕ_ d,d̅|^2 and |ϕ_-1|^2 |ϕ_ d,d̅|^2 (16 diagrams: 8 for +1 state  and 8 for -1 state). All these contributions have generally the same form so we write down only one of them, e.g. |ϕ_+1|^2 |ϕ_ d|^2 (there are 4 loops corresponding to this process):

    2/4T ∑_ k_1… k_4∑_Ω_1…Ω_3Δ̅_+1( k_2, k_1, Ω_1)Δ_+1( k_2,  k_3, Ω_2) Δ̅_ d( k_4,  k_3, Ω_3)Δ_ d( k_4,  k_1, Ω_1 + Ω_3 -Ω_2) 
    ×[-i Ω_1 + ε_c↓( k_2) - ε_v↑( k_1)  -i Ω_3 + ε_c↑( k_4) - ε_v↑( k_3) -i Ω_2 + ε_c↓( k_2) - ε_v↑( k_3)  - i Ω_1 +i Ω_2 -i Ω_3 + ε_c↑( k_4) - ε_v↑( k_1)].

In principle, one should also consider the three other expressions which have the same shape.

Furthermore, the spin-flip processes arise, i.e. ϕ_+1ϕ_-1ϕ̅_ dϕ̅_d̅ and ϕ_ dϕ_d̅ϕ̅_+1ϕ̅_-1 (there are 4+4 corresponding loops). For example the contribution from ϕ_ dϕ_d̅ϕ̅_+1ϕ̅_-1 in the form symmetrized regarding to the dispersion laws, is as follows:

    2/4T ∑_ k_1… k_4∑_Ω_1…Ω_3Δ̅_+1( k_2, k_1, Ω_1)Δ_d̅( k_2,  k_3, Ω_2) Δ̅_-1( k_4,  k_3, Ω_3)Δ_ d( k_4,  k_1, Ω_1 + Ω_3 -Ω_2) 
    ×[-i Ω_1 + ε_c↓( k_2) - ε_v↑( k_1)  -i Ω_3 + ε_c↑( k_4) - ε_v↓( k_3) -i Ω_2 + ε_c↓( k_2) - ε_v↓( k_3)  - i Ω_1 +i Ω_2 -i Ω_3 + ε_c↑( k_4) - ε_v↑( k_1)].


Now we can rewrite the fourth-order contributions in the c.m. variables, taking into account the saddle-point equations, and obtain the exciton-exciton interaction term and the saturation term with corrections due to dark excitons. Even though in the general case the 1s–exciton wavefunctions are different for Δ_± 1,2, for simplicity we will assume the band splitting negligible, so the hydrogen-like 1s–wavefunctions and the effective masses M_σ
will be the same. The effective exciton-photon action becomes 𝒮_ eff = 𝒮_0 + 𝒮_ int^b +𝒮_ int^d + 𝒮_ int^b-d + 𝒮_ int^ SF, where

    𝒮_0 =  ∫d  k/(2π)^2∫d Ω/2 π{∑_σ = ± 1,
     d,d̅C_σ( k, Ω)(-iΩ + E_ g + ħ^2  k^2/2 M_σ + E_ 1s^σ -μ)  C_σ( k, Ω) . 
    
         . + ∑_σ =± 1(Ψ_ ph, σ( k,Ω) [- i Ω + E_ ph( k) -μ] Ψ_ ph, σ( k,Ω) + ħΩ_ R/2[C_σ( k, Ω)Ψ_ ph, σ( k, Ω) +Ψ_ ph, σ( k, Ω) C_σ( k, Ω) ]) }

is the “bare” exciton-photon action without interactions,

    𝒮_ int^d + 𝒮_ int^b =  ∫d l_1d l_2d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3{∑_σ = ± 1,
     d,d̅g_ ex/2 C_σ( l_1, Ω_1) C_σ( l_2, Ω_2)  C_σ( l_3, Ω_3) C_σ( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3) .
    
        . -  g_ sat/2∑_σ =± 1[ C_σ( l_1, Ω_1) C_σ( l_2, Ω_2)  C_σ( l_3, Ω_3)Ψ_ ph,σ( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3 ) + c.c.] }

describe interactions between dark excitons and between bright excitons (within the same species), including bright-exciton-assisted exciton-photon coupling,

    𝒮_ int^b-d= ∫d l_1d l_2d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3∑_σ = ± 1{ g_ ex[C_ d( l_1, Ω_1) C_ d( l_2, Ω_2) + C_d̅( l_1, Ω_1) C_d̅( l_2, Ω_2)] . 
    ×  C_σ( l_3, Ω_3) C_σ( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3)  
    .  - g_ sat/2[C_ d( l_1, Ω_1) C_ d( l_2, Ω_2) + C_d̅( l_1, Ω_1) C_d̅( l_2, Ω_2)][C_σ( l_3, Ω_3)Ψ_ ph,σ( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3 ) +  c.c.]}

describes interactions between the dark and bright excitons of the same spin, including dark-exciton-assisted exciton-photon coupling and, finally,

    𝒮_ int^ SF = ∫d l_1d l_2d l_3/(2π)^6∫dΩ_1dΩ_2dΩ_3/(2π)^3{ g_ ex^ SF C_+1( l_1, Ω_1) C_d̅( l_2, Ω_2)  C_-1( l_3, Ω_3) C_ d( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3) .
        
    
        - g_ sat/2  C_d̅( l_2, Ω_2) C_ d( l_1-  l_2+ l_3, Ω_1 -Ω_2 +Ω_3) [C_+1( l_1, Ω_1) Ψ_ ph,-1( l_3, Ω_3) + Ψ_ ph, +1( l_1, Ω_1) C_-1( l_3, Ω_3)] 
     +  g_ ex^ SF C_ d( l_1, Ω_1) C_-1( l_2, Ω_2)  C_d̅( l_3, Ω_3) C_+1( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3)
        
    
        -  g_ sat/2 C_ d( l_1, Ω_1) C_d̅( l_3, Ω_3) [Ψ_ ph, -1( l_2, Ω_2) C_+1( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3) . +  C_-1( l_2, Ω_2) Ψ_ ph, +1( l_1 - l_2 + l_3, Ω_1 -Ω_2 +Ω_3)] }

corresponds to the spin-flip processes including also the dark-exciton-assisted decay.

It is important to note that other channels of pairing can also be taken into account in the same manner as done in Appendix <ref> where it was shown that only the off-diagonal density-like fields of the form Φ_c(v), σ,σ^'(k_1,k_2)= Ψ_c(v), σ(k_1) Ψ_c(v), σ^'(k_2)  contribute to the resulting exciton-exciton interaction (with the inclusion of spins there will be 8 off-diagonal density fields, 4 for conductance and 4 for the valence band). At the same time, the treatment of the 8 remaining diagonal density-like fields alters the Bethe-Salpeter equation (here, for different exciton fields, which is clearly seen from writing down the third-order expansion terms). When taking into account the additional off-diagonal fields, one needs to calculate 48 vertices of the same type: 24 of those lead to the interaction renormalization of the exciton fields within the same species, the other 24 diagrams screen the interaction between the dark and bright states. Importantly, the spin-flip process is not influenced by the off-diagonal density-like fields.
Thus in all the expressions obtained above for the effective action in the considered simplest case of zero temperature, rigid 1s–excitons and negligible band splitting, all the bright-bright, dark-dark, and bright-dark interaction constants are renormalized as derived in Eq. (<ref>). 
For spin-flip interactions, however, g_ ex remains unrenormalized, i.e.

    g_ ex^ SF = 2∫d p/(2π)^2∫d q/(2π)^2 V( q)χ( p)χ^*( p)χ^*( p - q)χ( p ).


To address corrections to the Rabi-splitting due interactions, we assume for simplicity that photons have polarization +1, and follow the procedure as described in Appendix <ref>, focusing on the +1 excitons:

    . δ𝒮_ eff/δC_+1|_C_+1=Ψ_ ex  
    Ψ_ ph, +1= Ψ_ ph
      C_-1      
    
         C_ d,d̅      
    Ψ_ ph, -1 = 0   = (E_ g+ E_1s^(+1)-μ) Ψ_ ex + ħΩ_ R/2Ψ_ ph +g_ ex |Ψ_ ex|^2 Ψ_ ex -g_ sat |Ψ_ ex|^2 Ψ_ ph - g_ sat/2Ψ_ ex^2 Ψ_ ph+ g_ ex (C_d̅C_d̅ + C_ dC_ d) Ψ_ ex
     
          - g_ sat/2(C_d̅C_d̅ + C_ dC_ d)Ψ_ ph +
         g_ ex^ SF C_-1C_ dC_d̅ -g_ sat/2 C_ dC_d̅Ψ_ ph,-1=0,

where the last term is equal zero since only the +1 polarization of photons is populated. Analogously,

    .δ𝒮_ eff/δΨ_ ph,+1|_
         C_+1=Ψ_ ex  
    Ψ_ ph, +1 = Ψ_ ph
      C_-1      
    
         C_ d,d̅      
    Ψ_ ph, -1 = 0   = (E_ ph^0-μ ) Ψ_ ph +ħΩ_ R/2Ψ_ ex - g_ sat/2 |Ψ_ ex|^2 Ψ_ ex - g_ sat/2(C_d̅C_d̅ + C_ dC_ d)Ψ_ ex - g_ sat/2 C_ dC_d̅C_-1=0 .

The (-1) exciton branch is assumed to be weakly populated compared to Ψ_ ex, so the terms containing C_-1 are neglected for clarity. We note that the last term describes the effective decay.

From the above equations, introducing the dark exciton density as n_ d = C_d̅ C_d̅+ C_ d C_ d, we express the bottoms of the UP and LP dispersion laws similarly as for the spinless case:

    μ_ LP(UP)(n_ ex^0, n_ d) = E_ g +E_1s^(+1) + g_ ex (n_ ex^0+ n_ d) +  Δ(n_ ex^0, n_d)/2±1/2√(Δ(n_ ex^0, n_ d)^2 + [ħΩ_ R(n_ ex^0, n_ d)]^2),

where Δ(n_ ex^0, n_ d) =Δ - g_ ex(n_ ex^0 + n_ d) and ħΩ_ R(n_ ex^0, n_ d) = √([ħΩ_ R - g_ sat(n_ ex^0+ n_ d)][ħΩ_ R - 3 g_ sat n_ ex^0- g_ sat n_ d]). The splitting between the polariton branches at the anticrossing point is

    min[E_ UP( k)- E_ LP( k)] = min√(Δ_ k^2(n_ ex^0, n_ d) + [ħΩ_ R(n_ ex^0, n_ d)]^2),

with Δ_ k(n_ ex^0, n_ d) = Δ - g_ ex(n_ ex^0 +n_ d)+ ħ^2  k^2/2 m_ ph - ħ^2  k^2/2 M. As one can see, the presence of dark exciton states effectively reduces the detuning and the Rabi-splitting with the increase of the dark-exciton density. 


00
NatRevQuantum A. Kavokin, T. C. H. Liew, C. Schneider, P. G. Lagoudakis, S. Klembt, and S. Hoefling, Polariton condensates for classical and quantum computing, https://doi.org/10.1038/s42254-022-00447-1Nat Rev Phys 4, 435 (2022) 

QFL I. Carusotto and C. Ciuti, Quantum fluids of light, https://link.aps.org/doi/10.1103/RevModPhys.85.299Rev. Mod. Phys. 85, 299 (2013)

estrecho E. Estrecho, T. Gao, N. Bobrovska, D. Comber-Todd, M. D. Fraser, M. Steger, K. West, L. N. Pfeiffer, J. Levinsen, M. M. Parish, T. C. H. Liew, M. Matuszewski, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya, Direct measurement of polariton-polariton interaction strength in the Thomas-Fermi regime of exciton-polariton condensation, https://link.aps.org/doi/10.1103/PhysRevB.100.035306Phys. Rev. B 100, 035306 (2019)

richard I. Frérot, A. Vashisht, M. Morassi, A. Lemaître, S. Ravets, J. Bloch, A. Minguzzi, and M. Richard, Bogoliubov excitations driven by thermal lattice phonons in a quantum fluid of light, https://doi.org/10.1103/PhysRevX.13.041058Phys. Rev. X 13, 041058 (2023)

menon  J. Gu, V. Walther, L. Waldecker, D. Rhodes, A. Raja, J. C. Hone, T. F. Heinz, S. Kéna-Cohen, T. Pohl, and V. M. Menon, Enhanced nonlinear interaction of polaritons via excitonic Rydberg states in monolayer WSe_2, https://doi.org/10.1038/s41467-021-22537-xNat Commun 12, 2269 (2021)

zhao J. Zhao, A. Fieramosca, K. Dini, R. Bao, W. Du, R. Su, Y. Luo, W. Zhao, D. Sanvitto, T. C. H. Liew, and Q. Xiong, Exciton polariton interactions in Van der Waals superlattices at room temperature, https://doi.org/10.1038/s41467-023-36912-3Nat. Commun 14, 1512 (2023)

stepanov P. Stepanov, A. Vashisht, M. Klaas, N. Lundt, S. Tongay, M. Blei, S. Höfling, T. Volz, A. Minguzzi, J. Renard, C. Schneider, and M. Richard, Exciton-exciton interaction beyond the hydrogenic picture in a MoSe_2 monolayer in the strong light-matter coupling regime, https://doi.org/10.1103/PhysRevLett.126.167401Phys. Rev. Lett. 126, 167401 (2021)

deng_const L. Zhang, F. Wu, Sh. Hou, Zh. Zhang, Y.-H. Chou, K. Watanabe, T. Taniguchi, S. R. Forrest, and H. Deng, Van der Waals heterostructure polaritons with moiré-induced nonlinearity, https://doi.org/10.1038/s41586-021-03228-5Nature 591, 61–65 (2021)


BEC_polaritons J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and Le Si Dang, Bose–Einstein condensation of
exciton polaritons, https://doi.org/10.1038/nature05131Nature 443, 409 (2006)

SF_polaritons A. Amo, J. Lefrère, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdré, E. Giacobino, and A. Bramati, Superfluidity of polaritons in semiconductor microcavities, https://doi.org/10.1038/nphys1364Nature Phys 5, 805–810 (2009)

yamamoto F. Tassone and Y. Yamamoto, Exciton-exciton scattering dynamics in a semiconductor microcavity and stimulated scattering into polaritons, https://doi.org/10.1103/PhysRevB.59.10830Phys. Rev. B 59, 10830 (1999)

glazov M. M. Glazov, H. Ouerdane, L. Pilozzi, G. Malpuech, A. V. Kavokin, and A. D’Andrea, Polariton-polariton scattering in microcavities: A microscopic theory, https://doi.org/10.1103/PhysRevB.80.155306Phys. Rev. B 80, 155306 (2009)

combescot M. Combescot, O. Betbeder-Matibet, and F. Dubin, The many-body physics of composite bosons, https://doi.org/10.1016/j.physrep.2007.11.003Phys. Rep. 463, 215 (2008)

parish_ehph1 J. Levinsen, G. Li, and M. M. Parish, Microscopic description of exciton-polaritons in microcavities, https://doi.org/10.1103/PhysRevResearch.1.033120Phys. Rev. Research 1, 033120 (2019)

schwendimannG. Rochat, C. Ciuti, V. Savona, C. Piermarocchi, A. Quattropani, and P. Schwendimann, Excitonic Bloch equations for a two-dimensional system of interacting excitons, https://doi.org/10.1103/PhysRevB.61.13856 Phys. Rev. B 61, 13856 (2000)

binder_scat R. Takayama, N.H. Kwong, I. Rumyantsev, M. Kuwata-Gonokami, and R. Binder, T-matrix analysis of biexcitonic correlations in the nonlinear optical response of semiconductor quantum wells, https://doi.org/10.1140/epjb/e20020051Eur. Phys. J. B 25, 445 (2002)

kira M. Kira, F. Jahnke, S. W. Koch, J. D. Berger, D. V. Wick, T. R. Nelson, Jr., G. Khitrova, and H. M. Gibbs, Quantum Theory of Nonlinear Semiconductor Microcavity Luminescence Explaining “Boser” Experiments, https://doi.org/10.1103/PhysRevLett.79.5170Phys. Rev. Lett. 79, 5170 (1997)

khurgin1 D. S. Citrin and J. B. Khurgin, Microcavity effect on the electron-hole relative motion in semiconductor quantum wells, https://doi.org/10.1103/PhysRevB.68.205325Phys. Rev. B 68, 205325 (2003)



PI4 P. R. Eastham and P. B. Littlewood, Bose condensation of cavity polaritons beyond the linear regime: The thermal equilibrium of a model microcavity, https://doi.org/10.1103/PhysRevB.64.235101Phys. Rev. B 64, 235101 (2001)

PI1 J. Keeling, P. R. Eastham, M. H. Szymanska, and P. B. Littlewood, Polariton condensation with localized excitons and propagating photons, https://doi.org/10.1103/PhysRevLett.93.226403Phys. Rev. Lett. 93, 226403 (2004)

PI3 M. H. Szymańska, J. Keeling, and P. B. Littlewood, Nonequilibrium quantum condensation in an incoherently pumped dissipative system, https://doi.org/10.1103/PhysRevLett.96.230602Phys. Rev. Lett. 96, 230602 (2006)

PI5 F. M. Marchetti, J. Keeling, M. H. Szymańska, and P. B. Littlewood, Absorption, photoluminescence, and resonant Rayleigh scattering probes of condensed microcavity polaritons, https://doi.org/10.1103/PhysRevB.76.115326Phys. Rev. B 76, 115326 (2007)

PI2 M. H. Szymańska, J. Keeling, and P. B. Littlewood, Mean-field theory and fluctuation spectrum of a pumped decaying Bose-Fermi system across the quantum condensation transition, https://doi.org/10.1103/PhysRevB.75.195331Phys. Rev. B 75, 195331 (2007)


PI6 J. Keeling, P. R. Eastham, M. H. Szymanska, and P. B. Littlewood, BCS-BEC crossover in a system of microcavity polaritons, https://doi.org/10.1103/PhysRevB.72.115320Phys. Rev. B 72, 115320 (2007)


parish_ehph3 S. S. Kumar, M. M. Parish, and J. Levinsen, Microscopic theory of excitons bound by light, https://doi.org/10.1103/PhysRevB.106.205414Phys. Rev. B 106, 205414 (2022)

keldysh L.V. Keldysh, A.N. Kozlov, Collective Properties of Excitons in Semiconductors, JETP 54, 978 (1968)

kiselev V.S. Babichenko and M.N. Kiselev, Superconductivoty in systems with excitonic instability, Journal of Moscow Phys. Soc. 2, 311 (1992)


tokatly A. A. Gorbatsevich and I. V. Tokatly, Phase transitions in an exciton semiconductor and in a superfluid low-density Fermi liquid: Bose-liquid approach. Three and two-dimensional systems. General description at finite temperatures,  JETP 108, 1723 (1995)






berkelbach_prb88 T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Theory of neutral and charged excitons in monolayer transition metal dichalcogenides, https://link.aps.org/doi/10.1103/PhysRevB.88.045318Phys. Rev. B 88, 045318 (2013)

chernikov_prl113 A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, B. Aslan, D. R. Reichman, M. S. Hybertsen, and T. F. Heinz, Exciton Binding Energy and Nonhydrogenic Rydberg Series in Monolayer 
WS_2, https://link.aps.org/doi/10.1103/PhysRevLett.113.076802Phys. Rev. Lett. 113, 076802 (2014)

keldysh_int L. V. Keldysh, Coulomb interaction in thin semiconductor and semimetal films, JETP Lett. 29, 658 (1979)

parish_pol O. Bleu, G. Li, J. Levinsen, and M. M. Parish, Polariton interactions in microcavities with atomically thin semiconductor layers, https://doi.org/10.1103/PhysRevResearch.2.043185Phys. Rev. Research 2, 043185 (2020)

LozYudson Yu. E. Lozovik and V. I. Yudson, On the ground state of the two-dimensional non-ideal bose gas, https://doi.org/10.1016/0378-4371(78)90170-XPhysica A 93, 493 (1978)

Castin_notes Y. Castin, Simple theoretical tools for low dimensional Bose gases, https://doi.org/10.1051/jp4:2004116004J. Phys. IV 116, 89 (2004)

trion-polaritons R. P. A. Emmanuele, M. Sich, O. Kyriienko, V. Shahnazaryan, F. Withers, A. Catanzaro, P. M. Walker, F. A. Benimetskiy, M. S. Skolnick, A. I. Tartakovskii, I. A. Shelykh, and D. N. Krizhanovskii, Highly nonlinear trion-polaritons in a monolayer semiconductor, https://doi.org/10.1038/s41467-020-17340-zNat Commun 11, 3589 (2020)

keldysh_1 L.V. Keldysh and A.N. Kopaev, Possible instability of the semimetallic state toward Coulomb interaction, Soviet Physics Solid State 6, 2219 (1965)

marchetti F. M. Marchetti, B. D. Simons, and P. B. Littlewood, Condensation of cavity polaritons in a disordered environment, https://doi.org/10.1103/PhysRevB.70.155327Phys. Rev. B 70, 155327 (2004)

lozovik_PI A. Elistratov, Yu. Lozovik, Coupled exciton-photon Bose condensate in path integral formalism, https://doi.org/10.1103/PhysRevB.93.104530Phys. Rev. B 93, 104530 (2016)

khurgin J.B. Khurgin, Excitonic radius in the cavity polariton in the regime of very strong coupling, https://doi.org/10.1016/S0038-1098(00)00469-5Solid State Communications, 117, 307 (2001)

Shahnazaryan2017 V. Shahnazaryan, I. Iorsh, I. A. Shelykh, and O. Kyriienko, Exciton-exciton interaction in transition-metal dichalcogenide monolayers, https://doi.org/10.1103/PhysRevB.96.115409Phys. Rev. B 96, 115409 (2017)

opala A. Opala, M. Furman, M. Król, R. Mirek, K. Tyszka, B. Seredyński, W. Pacuski, J. Szczytko, M. Matuszewski, and B. Piętka, Natural exceptional points in the excitation spectrum of a light–matter system, https://doi.org/10.1364/OPTICA.497170Optica 10, 1111 (2023)

wouters_biex M. Wouters, Resonant polariton-polariton scattering in semiconductor microcavities, https://doi.org/10.1103/PhysRevB.76.045319Phys. Rev. B 76, 045319 (2007)

vladimirova M. Vladimirova, S. Cronenberger, D. Scalbert, K. V. Kavokin, A. Miard, A. Lemaître, J. Bloch, D. Solnyshkov, G. Malpuech, and A. V. Kavokin, Polariton-polariton interaction constants in microcavities, https://doi.org/10.1103/PhysRevB.82.075301Phys. Rev. B 82, 075301 (2010)

parish_biex K. Choo, O. Bleu, J. Levinsen, and M. M. Parish, Polaronic polariton quasiparticles in a dark excitonic medium, https://doi.org/10.48550/arXiv.2312.00985arXiv:2312.00985

shelykh_biex A. Kudlis, I. A. Aleksandrov, M. M. Glazov, and I. A. Shelykh, Theory of biexciton-polaritons in transition metal dichalcogenide monolayers, https://doi.org/10.48550/arXiv.2402.09110arXiv:2402.09110

HS R.L. Stratonovich, On a Method of Calculating Quantum Distribution Functions, Soviet Physics, Doklady, 2, 416 (1958); J. Hubbard, Calculation of Partition Functions, https://doi.org/10.1103/PhysRevLett.3.77Phys. Rev. Lett. 3, 77 (1959)



