\section{Master equation and self-consistent approach}
We start discussing the case of incoherent injection of photons. For simplicity, we first neglect saturation effects, which are discussed in detail in Sec.~\ref{AppendixC}. The evolution of the cavity mode, determined by the Heisenberg equation of motion is
\begin{gather}
i\frac{d\hat a}{d t}=-\Delta\omega_c\hat a+\frac{\Omega}{\sqrt{N}}\sum_{i}\hat x_i-i\frac{\gamma_c}{2}\hat a+F,
\label{EqPh}
\end{gather}
now, it is convenient to re-write the cavity mode as $\langle \hat a\rangle=\sqrt{N}\alpha.$ In terms of this parameter, the steady-state can be written as
\begin{gather}
\alpha=\frac{1}{\Delta_c}\left(\frac{F}{\sqrt{N}}+\frac{\Omega}{N}\sum_{i}\langle \hat x_i\rangle\right)=\frac{1}{\Delta_c}(f+\Omega\langle x\rangle),
\label{alpha}
\end{gather}
where $\Delta_c=\Delta\omega_c+i\gamma_c/2$  and $f=F/\sqrt{N}$.

We use the solution for the steady-state of the photon field to define an effective Hamiltonian for the excitons, 
\begin{gather}
\hat H_{\text{local}}=\hat H_{X}+\sum_i(f_\Omega\hat x_i^\dagger+f_\Omega^*\hat x_i),
\end{gather}
where 
\begin{gather}
f_\Omega=\frac{\Omega}{\Delta_c}(f+\Omega\langle x\rangle),
\label{fOmega}
\end{gather}
here, $f_\Omega$ can be regarded as an effective exciton driving term. For relevant experimental parameters, the kinetic energy is heavily suppressed, that is, $U_X\gg t$. Therefore, we treat  exciton tunnelling at the mean-field level, $$-\sum_{ij}t_{ij}\hat x_i^\dagger\hat x_j\rightarrow -J\sum_i (\hat x_i^\dagger \langle x\rangle+ \hat x_i \langle x^\dagger\rangle)$$ with $J=zt,$ being $z$ the coordination number and $t$ the hopping coefficient. Then, we re-define $f_X$ to include at the mean-field level the hopping ratio $t$  
 \begin{gather}
 f_X=\frac{\Omega\,f}{\Delta_c}+\left(\frac{\Omega^2}{\Delta_c}-J\right)\langle \hat x\rangle,
 \end{gather}
the effects of exciton tunneling are discussed in detail in section ~\ref{NN}.

The evolution of the density matrix for the excitons is governed by the effective master equation 
\begin{gather}
\frac{d\hat\rho}{dt}=-i[\hat H_{\text{local}},\hat\rho]+\mathcal D_x[\hat\rho]=\mathcal L_x[\hat\rho],
\label{eqSs}
\end{gather}
where we should recall that the Hamiltonian $\hat H_{\text{local}}$ depends itself on the expectation $\langle \hat x\rangle=\text{Tr}(\hat x\hat\rho).$ Therefore, we adapt self-consistently the Lindblad operator. The self-consistent approach, which iterates on $\langle \hat x\rangle$ can be summarised as follows. We numerically diagonalise the local Lindblad operator taking a given value of  $\langle \hat x\rangle^{(i)},$ where $i$ denotes the $i$-th iteration step, which is employed in  Eq.~\ref{fOmega} to define the effective exciton Hamiltonian. We obtain the steady-state $\hat\rho^{(i)}$ which we use to evaluate the expectation value of $\hat x$. This defines $\langle \hat x\rangle^{(i+1)}=\text{Tr}(\hat x\hat\rho^{(i)}).$ We then replace $\langle \hat x\rangle^{(i+1)}$ in Eq.~\ref{fOmega} and repeat the iterative scheme. Numerically, we diagonalise the local Lindblad operator for $N_L=225$ corresponding to a cut-off for the exciton's level of $N=15$ and 200 iteration cycles unless stated explicitly otherwise. 

For $\Omega_{\text{sat}}=0,$ our numerics can be benchmarked against the analytical solutions which can be extended from the case of an isolated cavity with a Kerr-like non-linearity~\cite{Drummond1980}. The density matrix extended from ~\cite{Drummond1980}, can be written as a self-consistent equation
\begin{gather}
\label{DMatrix}
\hat \rho_{n,m}=\frac{(-2)^{n+m}}{\sqrt{n!m!}}\left(\frac{f_X}{U_X}\right)^n\left(\frac{f_X^*}{U_X}\right)^m\frac{\Gamma\left(\frac{-2\Delta_X}{U_X}\right)\Gamma\left(\frac{-2\Delta_X^*}{U_X}\right)}{\Gamma\left(\frac{-2\Delta_X}{U_X}+n\right)\Gamma\left(\frac{-2\Delta_X^*}{U_X}+m\right)}\frac{_0\mathcal F^{2}\left(\frac{-2\Delta_X}{U_X}+n,\frac{-2\Delta_X^*}{U_X}+m,4\left|\frac{f_x}{U_X}\right|^2\right)}{_0\mathcal F^{2}\left(\frac{-2\Delta_X}{U_X},\frac{-2\Delta_X^*}{U_X},8\left|\frac{f_X}{U_X}\right|^2\right)},
\end{gather}
where $_0\mathcal F^{2}(a,b;z)$ is the hypergeometric function defined as
\begin{gather}
_0\mathcal F^{2}(a,b;z)=\sum_{n=0}^{\infty}\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+n)\Gamma(b+n)}\frac{z^n}{n!},
\end{gather}
here, the self-consistency stems from the dependance of $f_X$ on $\text{Tr}(\hat\rho\,\hat x).$  Before comparing these approaches, we shortly provide details on the extended Gross-Pitaeksvii equation.
\subsection{Extended Gross-Pitaevskii equation and master equation}
We now turn our attention into the comparison between the semi-classical description of the excitons, the numerical self-consistent scheme, the density matrix in Eq.~\ref{DMatrix} obtained self-consistenly, and the non-interacting case. First, we consider non-interacting excitons $U_X=0$, in this case, the steady-state acquires the simple form of
\begin{gather}
\alpha=\frac{f}{\Delta_c}\frac{1}{1-\frac{\Omega^2}{\Delta_x\Delta_c}},\\ 
\langle \hat x \rangle=\frac{\Omega}{\Delta_x}\alpha.
\label{idealSS}
\end{gather}
On the other hand, after some algebra and introducing $\phi=\langle \hat x \rangle,$ the  extended driven-dissipative Gross-Pitaevskii for the steady-state reads can be written as
\begin{gather}
\left(-\Delta_X-J+\frac{\Omega^2}{\Delta_c}+U_X|\phi|^2\right)\phi+\frac{\Omega\,f}{\Delta_c}=0.
\end{gather}
We can now compare the coherence amplitude $|\phi|$ predicted from the different approximations. In Fig.~\ref{FigA1} we plot the coherence amplitude, the red squares correspond to the self-consistent approach described for the Lindblad operator $\mathcal L,$  while the solid blue line gives the benchmark for the density matrix in  Eq.~\ref{DMatrix}. This excellent agreement remarks the equivalence between these two numerical approaches. The black circles correspond to the eGP equation which predicts a bi-stability. Finally, the solid green line gives the solution for non-interacting excitons, namely Eq.~\ref{idealSS}. As discussed in the main text, for low intensities, the coherence amplitude is governed by Eq.~\ref{idealSS}. As soon as the interactions become relevant, the coherence amplitude obtained from Eq.~\ref{DMatrix} completely deviates from the eGP prediction illustrating the strongly interacting nature of the excitons, which cannot be explained by the eGP equation.


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\columnwidth]{FigA1}
\end{center}
\caption{Comparison between different approaches for the same value of the parameters than Fig.~3 (main text). The solid green line depicts the steady state for non-interacting excitons, the black circles correspond to the eGP equation, the red squares and solid blue line illustrate the equivalence between the self-consistent diagonalisation of the Lindblad operator and the iterative approach for the reduced density matrix in Eq.~\ref{DMatrix} respectively.} 
\label{FigA1}
\end{figure} 

\section{Nearest neighbour tunnelling and Saturation effects}

\subsection{Nearest neighbour tunnelling effects} 
\label{NN}
The hybridization of the excitonic band is accompanied by a strong suppression of their kinetic energy, this leads to a small hopping coefficient $t\ll U_X$ which can be understood at the mean-field level. 

The corrections due to nearest-neighbour tunnelling (NN) are expected to be small for
\begin{gather}\left|\frac{\Omega^2}{\Delta_c}\right|\gg zt=J,
\end{gather}
that is, as long as the second-order mediated tunnelling coefficient remains larger than the NN hopping ratio, we expect the tunnelling effects to be governed by the light-mediated tunnelling.  

To provide quantitative estimates of the NN hopping we consider a hopping coefficient $J/U_X=0.3.$  In this case, the ratio between the tunnelling amplitudes is  $\left|\frac{\Omega^2}{\Delta_c}\right|/J\approx 4$, for the light-matter coupling exemplified in the main text. In Fig.~\ref{FigAB1} we show the effects of NN tunnelling for a negative photon detuning $\Delta\omega_c=-1$ (red squares) compared against vanishing tunneling $J/U_X=0$ (blue circles). Fig.~\ref{FigAB1} shows that in this regime, the main features are only slightly affected by the NN hopping, in particular, the tunnelling mediated bi-stability is slightly enhanced by additional terms contributing to the hopping.  


\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.6\columnwidth]{FigB1}
\end{center}
\caption{Hysteresis and lower branch for $J/U_X=0.3$ (red squares) and $J/U_X=0$ (blue circles), which explicitly shows that the hysteresis and sudden jump are robust in presence of direct NN exciton tunnelling.  } 
\label{FigAB1}
\end{figure} 

\label{AppendixB}
\subsection{Saturation effects}
\label{AppendixC}
Saturation effects are included via an anharmonic light-matter coupling, such that the light-matter coupling reads as
\begin{gather}
\hat H_{l-m}=\frac{1}{\sqrt{N}}\sum_i \hat a^\dagger (\Omega-\Omega_{\text{sat}}\hat n_i)\hat x_i+\hat x_i^\dagger (\Omega-\Omega_{\text{sat}}\hat n_i)\hat a=\\ \nonumber
=\frac{\Omega}{\sqrt{N}}\sum_i \hat a^\dagger \left(1-\frac{\hat n_i}{n_\text{sat}}\right)\hat x_i+\hat x_i^\dagger \left(1-\frac{\hat n_i}{n_\text{sat}}\right)\hat a,
\end{gather}
the anharmonic term can be understood in terms of the suppression of the light-matter coupling $\Omega\rightarrow \Omega-\Omega_{\text{sat}}\hat n$. Then, saturation effects are relevant for densities of the order   $\langle \hat n \rangle \approx n_\text{sat}=\Omega_{\text{sat}}/\Omega.$ On the other hand, for  $ \langle n \rangle \ll n_\text{sat}$ the reduction of the light-matter coupling can safely be ignored. 



For two-dimensional gases, the reduction of the light-matter coupling $\Omega_{\text{sat}}=\Omega/(n_\text{sat}A)$ given in terms of the saturation density $n_{\text{sat}}=7/(16\pi a_B^2)$  with $a_B$ the exciton radius and the macroscopic area $A$~\cite{Cuiti2000} is a reminder that the intrinsic nature of the excitons becomes relevant when the inter-particle distance between excitons is of the order of the Bohr exciton radius. The spatial localisation of the hybrid excitons to the moir\'e lattice sites modifies the saturation condition but follows the same physical argument. However, the macroscopic area $A$ should be replaced by the ``local area" of the lattice site $A_M$. The local area is obtained by assuming that the spatial localisation of the excitons can be described by a periodic potential that can be approximated as a parabolic trap around each lattice site. Thus, around the minima of the lattice site, we consider a potential of the form $\frac{m \omega_M r_i^2}{2}.$ Therefore, the site dynamics of the exciton is governed by an harmonic oscillator Hamiltonian $$\hat H_{ho}=\frac{\hat p^2}{2m}+\frac{m \omega_M^2 \hat r^2}{2}.$$ For the ground-state, the classical return point $r_s=\frac{1}{\sqrt{m\omega_M}}$ defines the characteristic extension of the exciton wave function. Here, $m$ is the exciton mass, and $\omega_M $ is proportional to the electron tunnelling amplitude, responsible for the moir\'e lattice, typically of the order of $\omega_M=5-20\text{meV}$~\cite{Zhang2021}. We estimate  $r_s\approx 1-5\text{nm},$ whereas the typical exciton radius is $a_B\approx1\text{nm}$. The saturation density remains inversely proportional to the square of the exciton Bohr radius $n_s\propto 1/a_B^2$. Thus, the reduced light-matter coupling is estimated by $\Omega_{\text{sat}}/\Omega\propto (a_B/r_s)^2 \in \{0.1,1\}$. The suppression of the light-matter coupling sensible to one exciton per site is consistent with the experimental observation of moir\'e polaritons~\cite{Zhang2021}. This illustrates  the stark contrast between two-dimensional polariton gases and moir\'e polaritons which exhibit additional non-linearities even at the level of one and a few excitons per site~\cite{Zhang2021}. Here we take units where $\hbar=1$ and take $m=0.8m_e$ being $m_e$ the bare electron mass.


The anharmonic light-matter coupling modifies both the steady-state of the cavity mode, and the effective Hamiltonian for the excitons. The steady-state of the cavity acquires an explicit dependence on $\langle \hat n\hat x\rangle$  and reads as
\begin{gather}
\alpha=\frac{1}{\Delta_c}(f+\Omega\langle \hat x\rangle-\Omega_{\text{sat}}\langle \hat n\hat x\rangle),
\end{gather}
thus in addition to the dependence on the coherence $\phi=\langle \hat x\rangle$, the Hamiltonian develops an extra term consequence of the saturation of the exciton-photon coupling. In general, we should note that 
\begin{gather}
\Omega\langle \hat x\rangle-\Omega_{\text{sat}}\langle \hat n\hat x\rangle\neq (\Omega-\langle n\rangle\Omega_{\text{sat}})\langle\hat x\rangle,
\label{wrong}
\end{gather} thus we explicitly evaluate $\langle \hat n\hat x\rangle$. The right-side of Eq.~\ref{wrong} gives, however, intuition to understand the physical meaning of the saturation term. 
The inclusion of the anharmonic light-matter coupling renders the approach based on the analytical expression for the density matrix in Eq.~\ref{DMatrix} no longer valid. Therefore, the density matrix is now only obtained by means of the self-consistent method for the Lindblad super-operator.

To complement our discuss from the main text we show in Fig.~\ref{FigAC1} the coherence amplitude and purity as a function of the pump intensity and exciton detuning for the same parameters than in Fig.~2 (e) (main text). We show that for the low-densities, the lobular pattern persists and follows the discussion in the main text, where the transition from negligible site occupation to a finite number of exciton per site is accompanied by a drop of the purity and a non-monotonous behaviour of the coherence through the crossover.  In Fig.~\ref{FigAC1} we employ 100 iteration cycles and for the hysteresis branch we drive the system from an initial $f/\Omega=6$ (from above).


\begin{figure}[h!]
\begin{center}
\includegraphics[width=.65\columnwidth]{FigAC0}
\end{center}
\caption{Steady-state properties for the lower-branch for large saturation effects $\Omega_{\text{sat}}=\Omega/4.$ (Left) Coherence amplitude. (Right) Purity. We employ the same parameters than in Fig.~2 (e) (main text). } 
\label{FigAC1}
\end{figure} 

For negative exciton detuning, a bi-stable phase persists even when saturation effects become relevant,  the bi-stability is illustrated by cross-sections of the exciton number for fixed exciton detuning and as a function of $f/U_X$ in Fig.~\ref{FigAC2} for  $\Omega_{\text{sat}}=\Omega/4$, the red squares correspond to the hysteresis branch, and the blue circles to a $f/\Omega$ tuned from below.  Here we illustrate the persistence of the bi-stability for $\Delta\omega_X/U_X=0.6$ (left) and  $\Delta\omega_X/U_X=1.25$ (right).
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.35\columnwidth]{FigAC2}
\includegraphics[width=0.35\columnwidth]{FigAC2b}
\end{center}
\caption{Bi-stability  for $\Delta\omega_X/U_X=0.6$ (left) and  $\Delta\omega_X/U_X=1.25$ (right) and saturation effects of the order of $\Omega_{\text{sat}}=\Omega/4$. (Red squares depict the hysteresis branch, while blue circles correspond to the lower branch. } 
\label{FigAC2}
\end{figure} 

Finally, in Fig.~\ref{FigAC3} we compare the exciton number for different values of $\Omega_{\text{sat}}$ for $\Delta\omega_X/U_X=1.25,$ the hysteresis branch.  The red squares correspond to $\Omega_{\text{sat}}/\Omega=0.1$, the blue circles $\Omega_{\text{sat}}/\Omega=0.25$, and the black triangles $\Omega_{\text{sat}}/\Omega=0.5$. In this case, even for $\Omega_{\text{sat}}/\Omega=0.25$ saturation effects become prominent, the abrupt jump for small   $\Omega_{\text{sat}}/\Omega$ is visibly smoothen, the exciton number decreases, and the jump is shifted to larger values of $f/\Omega$. 
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.35\columnwidth]{FigAC3}
\end{center}
\caption{Exciton number for several values of $\Omega_{\text{sat}}$.  The red squares correspond to $\Omega_{\text{sat}}/\Omega=0.1$, the blue circles to $\Omega_{\text{sat}}/\Omega=0.25$, and the black triangles to $\Omega_{\text{sat}}/\Omega=0.5$ } 
\label{FigAC3}
\end{figure} 
\section{Lasers: From single-photon to multi-photon resonances}
\subsection{Effective two-level system}
Let us turn our attention into the details on the single- and multi-photon resonances in the scenario with exciton gain. For simplicity, we start discussing the limit where only empty moir\'e lattice sites can be driven, in this case, the gain and losses processes are described by\begin{gather}
\mathcal D_x[\hat\rho]=\Gamma_x\mathcal D[\hat y^\dagger]\hat\rho+\gamma_x\mathcal D[\hat x]\hat \rho,
\end{gather}
here, the first term accounts for the driving of empty sites, such that $\hat y$ and $\hat y^\dagger$ are exciton operators restricted to $N_x=0,1$ subspace. Now, we solve the equation $\mathcal L[\hat\rho]=0$ in a rotating frame with frequency $\omega.$ We define $\Delta\omega_X=\omega-\omega_X$ and $\Delta\omega_c=\omega-\omega_c$, where the photon energy $\omega$ is also left to be obtained. In the limit of $U_X\gg \Omega,$ the system can be regarded as a two-level system with  equations for the steady-state
\begin{gather}
0=\gamma \rho_{11}+\frac{4 \gamma_c\rho_{01}\rho_{10} \Omega ^2}{\gamma_c^2+4 (\delta +\omega )^2}-\Gamma_x\rho_{00},\\
0=\frac{2\rho_{01} \Omega ^2 (\rho_{11}-\rho_{00})}{\gamma_c+2 i (\delta +\omega )}-\frac{1}{2}\rho_{01} (\gamma +\Gamma_x+2 i  \omega ),\\ 
0=\frac{2\rho_{10} \Omega ^2 (\rho_{11}-\rho_{00})}{\gamma_c-2 i (\delta + \omega )}-\frac{1}{2}\rho_{10} (\gamma +\Gamma_x-2 i  \omega ),
\end{gather}
with $\delta=\omega_X-\omega_c$. In this case, there is a trivial solution with $\rho_{00}=\gamma_x/(\gamma_x+\Gamma_x)$ and $\rho_{01}=\rho_{10}=0$. A non-trivial solution appears for $\Omega>\Omega_c$ and $\Gamma_x>\gamma_x$, being $\Omega_c=(\gamma_x+\Gamma_x)\sqrt{\frac{\gamma_c}{4(\Gamma_x-\gamma_x)}},$ with  
\begin{gather}
\rho_{11} =\frac{1}{2}\left(1+\frac{\gamma_c(\gamma_x+\Gamma_x)}{4\Omega^2}\right)+\frac{1}{2}\left(\delta^2\frac{\gamma_c(\gamma_x+\Gamma_x)}{2(\gamma_x+\Gamma_x+\gamma_c)\Omega^2}\right)\\  \label{EnT}
\rho_{01}=\frac{i \sqrt{\gamma_c} \sqrt{(\gamma +\gamma_c+\Gamma_x)^2+4 \delta ^2} \sqrt{\gamma_c (\gamma +\Gamma_x)^2 \left((\gamma +\gamma_c+\Gamma_x)^2+4 \delta ^2\right)+4 \Omega ^2 (\gamma -\Gamma_x) (\gamma +\gamma_c+\Gamma_x)^2}}{4 \sqrt{2} \sqrt{\Omega ^4 (\gamma +\gamma_c+\Gamma_x)^4}},\\ \label{phEr}
\omega=\omega_X-\frac{\gamma_x+\Gamma_x}{\gamma_x+\Gamma_x+\gamma_c}\delta.
\end{gather}

In Fig.~\ref{FigAC4} we show the expectation value of $\hat x$ as a function of $\Omega$ for the two-level system in presence of exciton driving. The red dots correspond to the numerical implementation of the self-consistent approach. The solid blue line illustrates the analytical solution in Eq.~\ref{EnT}, while the dashed black line gives the onset for the non-trivial state. The onset $\Omega_c$ increases with $\gamma_x\rightarrow\Gamma_x$.  In this case, we find the photon energy self-consistently by solving the equivalent of Eq.~\ref{EqPh}. In general the photon energy slightly deviates from $\omega_c$ as shown by Eq.~\ref{phEr} depending on the cavity detuning $\delta$, and the damping and driving rates.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.6\columnwidth]{Fig2L}
\end{center}
\caption{Expectation value $\langle \hat x\rangle$ as a function of $\Omega$ for $\Gamma_x/\gamma_x=6,$ and $\Gamma_x=\gamma_c$. The red dots correspond to the self-consistent exact diagonalization scheme, the blue line the analytical solution in Eq.~\ref{EnT}. The vertical black dashed line gives $\Omega_c$. } 
\label{FigAC4}
\end{figure} 

For the two-photon resonance, we cannot longer derive an analytical expression for the photon amplitude, instead, we solve for the Lindblad operator for the exciton field $\mathcal L=\mathcal L(\omega,\alpha),$ where the photon energy $\omega$ and amplitude $\alpha$ are variables obtained self-consistently.  We stress that this self-consistent approach recovers the exact solution for the two-level system in Fig.~\ref{FigAC4} as well as, for the coherent injection of photons in absence of saturation effects, as discussed previously.  


\subsection{Three-level effective Hamiltonian: Rate equations}
To understand better the two-photon resonance, we derive the rate equations for an effective three-level system. We start by solving for the off-diagonal elements of the exciton's density matrix, which evolution is given by
\begin{gather}
\frac{d \rho_{01}}{dt}=-i \Omega  \left(\sqrt{2} \alpha \rho_{02}+\alpha^{*} (\rho_{00}-\rho_{11})\right)+\sqrt{2} \gamma_x \rho_{12}+\frac{1}{2} i\rho_{01} (U_X+i (\gamma_x +3 \Gamma_x),\\
\frac{d \rho_{02}}{dt}=i \alpha^{*} \Omega  \left(\rho _{12}-\sqrt{2}\rho_{01}\right)-\gamma_x \rho_{02}-\frac{\Gamma_x\rho_{02}}{2},\\ 
\frac{d \rho_{12}}{dt}=i \Omega  \left(\alpha \rho_{02}+\sqrt{2} \alpha^{*} (\rho_{22}-\rho_{11})\right)-\frac{3 \gamma_x \rho_{12}}{2}+\Gamma_x \left(\sqrt{2}\rho_{01}-\rho_{12}\right)-\frac{i\rho_{12} U}{2}.
\end{gather}
for the steady-state, and after some algebra, we find  a solution for the off-diagonal elements of the density matrix in terms of only its diagonal elements. These solutions are
\begin{gather}
\rho_{01}=\frac{2 \alpha^{*} \Omega  \left(-4 n_\alpha \Omega ^2 (\rho_{00}+\rho_{11}-2\rho_{22})+(2 \gamma_x -\Gamma_x (\gamma  (3\rho_{00}+\rho_{11}-4\rho_{22})+2 \Gamma_x (\rho_{00}-\rho_{11})+i U_X (\rho_{00}-\rho_{11}))\right)}{i (2 \gamma_x -\Gamma_x \left(3 \gamma_x ^2+3 \gamma_x  \Gamma_x+6 \Gamma_x^2+U_X^2-2 i \gamma  U+i \Gamma_x U_X\right)+4 \alpha  \alpha^{*} \Omega ^2 (U_X-3 i (\gamma_x +\Gamma_x)},\\ 
\rho_{12}=\frac{2 i \sqrt{2} \alpha^{*} \Omega  \left(-4 n_\alpha \Omega ^2 (\rho_{00}+\rho_{11}-2\rho_{22})+(2 \gamma_x -\Gamma_x (\Gamma_x (2\rho_{00}+\rho_{11}-3\rho_{22})+(\rho_{11}-\rho_{22}) (\gamma_x -i U_X))\right)}{4 \alpha  \alpha^{*} \Omega ^2 (3 (\gamma_x +\Gamma_X+i U_X)-(2 \gamma_x-\Gamma_x \left(3 \gamma_X ^2+3 \gamma_x  \Gamma_x+6 \Gamma_x^2+U_X^2-2 i \gamma_x  U_X+i \Gamma_x U_X\right)},\\
\rho_{02}=-\frac{4 \sqrt{2} (\alpha^{*})^2 \Omega ^2 (U_X(\rho_{00}-\rho_{22})-3 i (\gamma_x  (\rho_{00}-\rho_{22})+\Gamma_x (\rho_{22}-\rho_{11})))}{i (2 \gamma_x -\Gamma_x \left(3 \gamma_x ^2+3 \gamma_x  \Gamma_x+6 \Gamma_x^2+U^2-2 i \gamma_x U+i \Gamma_x U_X\right)+4 \alpha  \alpha^{*} \Omega ^2 (U_X-3 i (\gamma_x +\Gamma_x)}.
\end{gather}
To make progress, we make some simplifications, first we work in the regime where $\Gamma_x\gg \gamma_x$ and fix $\gamma_x=0,$ then we expand the off-diagonal elements in powers of  $\tilde\Omega=\Omega/U_X$ and retain terms up to $\tilde \Omega^3$. Thus, we obtain
\begin{gather}
\rho_{01}=2 \alpha^{*} \tilde\Omega  \left(\frac{8 n_\alpha (2 \tilde\Gamma_x+i) \tilde\Omega ^2 (-3 i \tilde\Gamma_x (\rho_{11}-\rho_{22})-\rho_{00}+\rho_{22})}{\tilde\Gamma_x (1+\tilde\Gamma_x (6 \tilde\Gamma_x+i))^2}-\frac{i (\rho_{00}-\rho_{11})}{3 \tilde\Gamma_x-i}\right),\\
\rho_{12}=\frac{8 \sqrt{2} n_\alpha\,\alpha^{*} (\tilde\Gamma_x+i) \tilde\Omega ^3 (-3 i \tilde\Gamma_x (\rho_{11}-\rho_{22})-\rho_{00}+\rho_{22})}{\tilde\Gamma_x (1+\tilde\Gamma_x (6 \tilde\Gamma_x+i))^2}+\frac{2 \sqrt{2} \alpha^{*} \tilde\Omega  (-i \tilde\Gamma_x (2\rho_{00}+\rho_{11}-3\rho_{22})-\rho_{11}+\rho_{22})}{1+\tilde\Gamma_x (6 \tilde\Gamma_x+i)},\\
\rho_{02}=\frac{4 \sqrt{2} (\alpha^{*})^2 \tilde\Omega ^2 (3 \tilde\Gamma_x (\rho_{11}-\rho_{22})-i\rho_{00}+i\rho_{22})}{\tilde\Gamma_x (1+\tilde\Gamma_x (6 \tilde\Gamma_x+i))},
\end{gather}
where we introduced $\tilde\Gamma_x=\Gamma_x/U_X.$ For the two-photon process, the single-particle transitions are far detuned, and for two-photon gain purposes, the linear terms in $\tilde \Omega$ can be neglected. One can further simplify these expressions by taking $\tilde \Gamma_x\ll 1$. Then, the elements for the off-diagonal density matrix take the following form
\begin{gather}
\label{g2rho}
\rho_{01}\approx \frac{16 i n_\alpha  \alpha^{*} \tilde\Omega ^3 (\rho_{22}-\rho_{00})}{\tilde\Gamma_x}+16 \alpha^*n_\alpha \tilde\Omega ^3 (-4 \rho_{00}+3 \rho_{11}+\rho_{22})\approx\frac{16 i n_\alpha  \alpha^{*} \Omega ^3 (\rho_{22}-\rho_{00})}{U_X^2\Gamma_x},\\ \label{g2rho1}
\rho_{12}\approx \frac{8 i \sqrt{2} n_\alpha  \alpha^{*}\tilde\Omega ^3 (\rho_{22}-\rho_{00})}{\tilde\Gamma_x}+24 \sqrt{2} n_\alpha\alpha^* \tilde\Omega ^3 (\rho_{11}-\rho _{00})\approx\frac{8 i \sqrt{2} n_\alpha  \alpha^{*}\Omega ^3 (\rho_{22}-\rho_{00})}{U_X^2\Gamma_x},\\ \label{g2rho2}
\rho_{02}\approx-\frac{4 i \sqrt{2} (\alpha^{*})^2 \tilde\Omega ^2 (\rho_{00}-\rho_{22})}{\tilde\Gamma_x}-4 \left(\sqrt{2} (\alpha^*)^2 \tilde\Omega ^2 (\rho_{00}-3\rho _{11}+2 \rho_{22})\right)\approx-\frac{4 i \sqrt{2} (\alpha^{*})^2 \Omega ^2 (\rho_{00}-\rho_{22})}{U_X\Gamma_x},
\end{gather}
where we have made the natural assumption that $\Omega\gg \Gamma_x$.

The equations for the diagonal elements of the density matrix are given by
\begin{gather}
\frac{d \rho_{00}}{dt}=-i \Omega  (\alpha \rho_{01}-\alpha^{*} \rho_{10})+\gamma_x \rho_{11}-\Gamma_x\rho_{00},\\
\frac{d \rho_{11}}{dt}=i \Omega  \left(\alpha(\rho_{01}-\sqrt{2}\rho_{12})+\alpha^{*}(- \rho_{10}+\sqrt{2} \rho_{21})\right)-\gamma_x  (\rho_{11}-2\rho_{22})+\Gamma_x (\rho_{00}-2\rho_{11}),\\
\frac{d\rho_{22}}{dt}=i \sqrt{2} \Omega  (\alpha \rho_{12}-\alpha^{*} \rho_{21})-2 \gamma_x \rho_{22}+2 \Gamma_x\rho_{11},
\end{gather}
we then replace the simplified expressions for the off-diagonal elements found in Eq.~\ref{g2rho}-\ref{g2rho2} to write the rate equations as follows 
\begin{gather}
\frac{d \rho_{00}}{dt}=+32\frac{ n^2_\alpha  \Omega ^4 (\rho_{22}-\rho_{00})}{U_X^2\Gamma_x}+\gamma_x \rho_{11}-\Gamma_x\rho_{00}+=8\frac{g_2^2}{\Gamma_x} n^2_\alpha (\rho_{22}-\rho_{00})+\gamma_x \rho_{11}-\Gamma_x\rho_{00},\\
\frac{d \rho_{11}}{dt}=-\gamma_x  (\rho_{11}-2\rho_{22})+\Gamma_x (\rho_{00}-2\rho_{11}),\\
\frac{d\rho_{22}}{dt}=-32\frac{ n^2_\alpha  \Omega ^4 (\rho_{22}-\rho_{00})}{U_X^2\Gamma_x}-2 \gamma_x \rho_{22}+2 \Gamma_x\rho_{11}=-8\frac{g_2^2}{\Gamma_x}n^2_\alpha  (\rho_{22}-\rho_{00})-2 \gamma_x \rho_{22}+2 \Gamma_x\rho_{11}.
\end{gather}
For $\gamma_x=0,$ the equation for $\rho_{11}$ simply gives $\rho_{00}=2\rho_{11}$ and yields the following equations
\begin{gather}
\frac{d \rho_{00}}{dt}=-\Gamma_x\rho_{00}+8\frac{g_2^2}{\Gamma_x} n^2_\alpha (\rho_{22}-\rho_{00}),\\
\frac{d\rho_{22}}{dt}= \Gamma_X\rho_{00}-8\frac{g_2^2}{\Gamma_x} n^2_\alpha (\rho_{22}-\rho_{00}),
\end{gather}
it is straightforward to note that  $\frac{d \rho_{00}}{dt}=-\frac{d \rho_{22}}{dt}.$ 

Finally, for the photon number we have the following equation
\begin{gather}
\frac{d n_\alpha}{dt}=-\gamma_c n_\alpha+i\Omega[\alpha^*(\rho_{10}+\sqrt{2}\rho_{21})-\alpha(\rho_{01}+\sqrt{2}\rho_{12})]\approx \\ \nonumber
\approx -\gamma_c n_\alpha+16\frac{g_2^2}{\Gamma_x} n^2_\alpha (\rho_{22}-\rho_{00}),
\end{gather}
from these equations  the two-photon emission rate is ready to be identified as   $A_{2}=8g_2^2/\Gamma_X.$ The equations can be written in its simplest form of
\begin{gather}
\label{ratee}
\frac{d \rho_{22}}{dt}\approx\Gamma_x\rho_{00}-A_2n_\alpha^2(\rho_{22}-\rho_{00}), \\ \nonumber
\frac{d n_\alpha}{dt}\approx-\gamma_cn_\alpha+2A_2 n_\alpha^2(\rho_{22}-\rho_{00}).
\end{gather}
 

The rate equations in Eq.~\ref{ratee} accept always a trivial solution with vanishing photon number $n_\alpha=0,$ and inverted population distributions $\rho_{00}=\rho_{11}=0$ and $\rho_{22}=1.$ However, beyond a critical two-photon emission rate a non-trivial steady state arises, this  critical $A_2^c$ is given by
\begin{gather}
A_2^{c}=\frac{5\gamma_c^2}{2\Gamma_x}, 
\end{gather}
and the photon number at the onset of two-photon gain is independent of $A_2$ 
\begin{gather}
n_\alpha=\frac{2\Gamma_x}{5\gamma_c}.
\end{gather}

One can further re-write the threshold for two-photon gain in terms of the two-photon coupling 
\begin{gather}
\frac{g^c_2}{\gamma_c}\approx \sqrt{\frac{5}{16}}\approx 0.55,
\end{gather}
where $g^c_2$ denotes the critical two-photon coupling.

In general, for finite $\gamma_x$, the trivial state has populations  given by  $\rho_{00}=\gamma_x^2/(\gamma_x^2+\gamma_x\Gamma_x+\Gamma_x^2)$ and $\rho_{11}=\gamma_x\Gamma_x/(\gamma_x^2+\gamma_x\Gamma_x+\Gamma_x^2).$  On the other hand, the appearance of the non-trivial state is slightly modified to
\begin{gather}
A_2^{c}=\frac{\gamma_c^2(5\Gamma_x+4\gamma_x)(\Gamma_x^2+\Gamma_x\gamma_x+\gamma_x^2)}{2(\Gamma_x^2-\gamma_x^2)^2},
\end{gather}
however, the photon number remains independent of $A_2$,
\begin{gather}
n_\alpha=2\frac{\Gamma_x^2-\gamma_x^2}{\gamma_c(5\Gamma_x+4\gamma_x)},
\end{gather}
note that solutions only appear for $\Gamma_X>\gamma_X$ as physically expected.

Finally, we shall remember that our numerical solutions work in the regime where $U_X\propto \Omega,$ that is, the exciton interactions  are of the order of magnitude of the light-matter coupling, leading to differences between the rate equations and the full numerical solutions. In addition, we have assume that the photon energy matches the cavity energy, $\omega=\omega_c$, contrary to the full numerics where $\omega$ is allowed to adapt self-consistently. However, the rate equations provide a valuable guide to understand the two-photon gain.
\end{widetext}
\end{document}
  
 


 
 
 
 
 

 
 
 


