../imgs/

<EMAIL_ADDRESS>
# Dynamics and condensation of polaritons in an optical nanocavity coupled to
two-dimensional materials

Maria Vittoria Gurrieri  [  Department of Photonics Engineering, Technical
University of Denmark, 2800 Kgs. Lyngby, Denmark  NanoPhoton - Center for
Nanophotonics, Technical University of Denmark, Ørsteds Plads 345A, DK-2800
Kgs. Lyngby, Denmark  Emil V. Denning  Department of Photonics Engineering,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  NanoPhoton -
Center for Nanophotonics, Technical University of Denmark, Ørsteds Plads 345A,
DK-2800 Kgs. Lyngby, Denmark  Nichtlineare Optik und Quantenelektronik,
Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin,
Germany  Kristian Seegert  Department of Photonics Engineering, Technical
University of Denmark, 2800 Kgs. Lyngby, Denmark  NanoPhoton - Center for
Nanophotonics, Technical University of Denmark, Ørsteds Plads 345A, DK-2800
Kgs. Lyngby, Denmark  Philip T. Kristensen  Department of Photonics
Engineering, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark
NanoPhoton - Center for Nanophotonics, Technical University of Denmark,
Ørsteds Plads 345A, DK-2800 Kgs. Lyngby, Denmark  Jesper Mørk  Department of
Photonics Engineering, Technical University of Denmark, 2800 Kgs. Lyngby,
Denmark  NanoPhoton - Center for Nanophotonics, Technical University of
Denmark, Ørsteds Plads 345A, DK-2800 Kgs. Lyngby, Denmark

###### Abstract

We present a comprehensive investigation of the light-matter interaction
dynamics in two-dimensional materials coupled with a spectrally isolated
cavity mode in the strong coupling regime. The interaction between light and
matter breaks the translational symmetry of excitons in the two-dimensional
lattice and results in the emergence of a localized polariton state. Employing
a novel approach involving transformation to exciton reaction coordinates, we
derive a Markovian master equation to describe the formation of a macroscopic
population in the localized polariton state. Our study shows that the
construction of a large-scale polariton population is affected by correction
terms addressing the breakdown of translational symmetry. Increasing the
spatial width of the cavity mode increases the Coulomb scattering rates while
the correction terms saturate and affect the system’s dynamics progressively
less. Tuning the lattice temperature can induce bistability and hysteresis
with different origins than those recognized for quantum wells in larger
microcavities. We identify a limit temperature $T_{\mathrm{l}}$ as a key
factor for stimulated emissions and forming a macroscopic population,
enriching our understanding of strong coupling dynamics in systems with
extreme confinement.

††preprint: APS/123-QED

##  I Introduction

In the study of light-matter interaction, strong coupling refers to the
condition where the coupling strength between photons and excitations in a
material dominates over dissipative processes [1]. This gives rise to many
fascinating phenomena because of the formation of polaritons, quasi-particles
that combine the properties of light and matter. In particular, when the
excitons in a quantum well or a two-dimensional (2D) semiconductor [2]
hybridize with the light field confined in an optical cavity, a condensate of
polaritons can be observed [3, 4, 5, 6, 7]. One mechanism behind the formation
of a macroscopic population of quasi-particles has been traced back to
stimulated scattering due to the quasi-bosonic nature of polaritons. Many
experiments have been conducted to study 2D and zero-dimensional (0D)
polariton systems [8], where these particles are confined in either one or
three dimensions. Theoretical studies in quantum wells describe polariton
dynamics through scattering processes conserving both energy and momentum [9,
10, 11].

Figure 1: Ground state polariton occupation as a function of pump power for
different in-plane confinement lengths, $L$, of the optical cavity. The curves
are computed by adapting the model from Ref. [10] to the case of a system with
a single optical mode. For each confinement length, the star indicates the
Mott transition.

These models predict that by increasing the confinement, i.e. reducing the in-
plane cavity length, the threshold to quantum degeneracy effects is reached at
progressively smaller pump powers. Such behaviour is illustrated in Fig. 1,
where the ground state polariton occupation number $N_{0}$ is plotted as a
function of pump power $p_{x}$ for different confinement lengths of a 2D
structure. The reduction in condensation threshold with in-plane confinement
motivates the interest in studying polariton condensation in systems with
extreme confinement. Current theoretical models, however, are derived for to
2D structures with transverse extent on the order of micrometres and larger.
Therefore, the usefulness and validity of these models are not guaranteed when
the optical cavity is scaled down to a transverse extent so small that the
system response is dominated by a single or a few resonances rather than a
continuum of modes. In particular, the assumption of translational symmetry
used in previous works [10] breaks down.

Figure 2: Schematic representation of a sheet of 2D material coupled to an
optical nanocavity, with a single optical mode localized at the centre of the
bowtie. Strong light-matter coupling between the resonant electromagnetic
field and the excitons in the 2D sheet leads to the formation of a localized
polariton state.

In this work, we study light-matter interaction in a 2D semiconductor with a
direct band gap placed on top of an optical nanocavity, as illustrated in Fig.
2. Within this system, the 2D semiconductor excitations form a particle
reservoir composed of thermalized excitons that can strongly couple to the
resonant electromagnetic field of the nanocavity. The nanocavity design is
taken from Ref. [12], and we note that a similar cavity was recently
experimentally demonstrated to implement deep sub-wavelength confinement [13].
Our goal is to reduce the threshold for detecting a macroscopic polariton
population, providing, at the same time, an extensive reservoir that can pump
particles into the condensate. The 2D structure coupled to the single cavity
mode is modelled as an open quantum system in the framework of exciton
reaction coordinates (ERC) [14, 15, 16], where the excited particles in the 2D
material behave as interacting bosons. The dynamics of the polariton
population are described by combining the formulation of a localized exciton
state in the ERC formalism with the rate-equation approach previously
developed by Tassone $et\ al.$ [9] and Porras $et\ al.$ [10]. The gain
mechanism, which is induced by stimulated scattering, is responsible for a
macroscopic population and the formation of a condensate in the lower
polariton state. This condensate, in turn, is responsible for coherent light
emissions and requires much lower excitation power than a standard laser. The
reduced threshold for coherent emissions is attributed to the fact that
condensation does not necessitate population inversion [17, 18].

The article is organized as follows. In Sec. II, the theoretical model is set
up. First, the ERC formalism is introduced. This enables a description in
terms of a system where a single exciton is coupled to one mode of the
electromagnetic resonator and a reservoir of residual excitons. The system’s
time evolution is described by a Markovian master equation from which a set of
rate equations for macroscopic observables is derived. In Sec. III, the main
results predicted by the theory are presented and explained. Specific
attention is given to the understanding of how the temperature, the in-plane
confinement, and the non-resonant pumping affect the dynamics. Sec. IV
presents the conclusions and further discussions.

##  II Theoretical framework

This section describes the mathematical tools employed to study the
interaction between light in a general electromagnetic resonator dominated by
a single resonance and excitons in an adjacent 2D system, as illustrated in
Fig. 2 for the case of an optical nanocavity. In a semiconductor, the Coulomb
interaction between electrons and holes results in the creation of excitons,
which can freely move in a translationally invariant medium. However, when
introducing the electromagnetic resonator, the interaction between light and
matter disrupts the translational invariance, giving rise to localized
excitonic states with a centre-of-mass wave function exactly matching the
electromagnetic field profile [19].

Figure 3: Schematic representation of the construction of the exciton reaction
coordinate formulation. (a) The localized resonant field, with $a$ the
operator describing the quantum properties of the resonant field, is coupled
to a continuum of exciton modes with momentum $\mathbf{k}$, described by the
annihilation operators ${b_{\mathbf{k}}}$. (b) The strength of the
light–matter interaction is described through the coupling strength $G_{0}$
between the resonant field and a single collective exciton reaction coordinate
with annihilation operator $B_{0}$, which, in turn, is coupled to an
environment of residual exciton modes described by the operators $B_{i}$. (c)
In the strong coupling regime, the Hopfield transformation [20, 21] enables a
description of the system through polariton quasi-particles with annihilation
operators $p_{\mathrm{lp/up}}$, forming the upper and lower polariton branch.
For the present applications, the former is energetically decoupled from the
dynamics and is therefore neglected.

If the light-matter coupling is sufficiently strong, it will give rise to a
spatially localized polariton state as illustrated in Fig. 3. The system is
initially described by coupling the resonant electromagnetic field to all
excitons with different in-plane momenta. In a change of basis, this is
equivalent to the case where the electromagnetic field is coupled to only a
single exciton reaction coordinate which, in turn, is coupled to a bath of
residual exciton states. In the strong coupling regime, finally, a Hopfield
transformation enables the identification of discrete lower and upper
polariton states. As shown below, the population of the lower polariton branch
can be modelled using a Born-Markov master equation, which considers the
interaction with a continuous reservoir of excitons. This, in turn, enables a
rate equation description for the dynamics of the lower polariton population.
A complete description of the dynamics requires two additional equations that
capture the evolution of the exciton density $n_{x}$, and the energy density
$e_{x}$ in the surrounding 2D material. The scattering processes between the
particles also modify the temperature and can lead to heating the surrounding
lattice.

###  A Exciton reaction coordinate and strong coupling

The model is based on Wannier-Mott exciton states and their interactions with
a resonant electromagnetic field described by a single quasi-normal mode [22,
23, 24, 25]. For the reference calculations and the simulation in Sec. IV, we
consider a monolayer of transition metal dichalcogenide (TMDC). These
materials present a direct band-gap at the $K$ and $K^{\prime}$ points of the
reciprocal lattice space [26, 27]. Their advantage, compared to semiconductor
quantum wells such as gallium arsenide, relates to the properties of monolayer
TMDC excitons [28, 29, 30]. Their large binding energy makes them stable and
enables the study of polariton condensation even close to room temperature
[31]. Also, their small Bohr radius provides access to higher exciton
densities before the Mott-transition is reached.

Wannier-Mott excitons are described via the relative position of the electron
and hole and their centre-of-mass, which is free to move throughout the entire
sheet of the 2D material. In practice, the centre of mass motion is described
by a plane wave of momentum $\mathbf{k}$. Being composed of electron-hole
pairs, excitons are represented by operators with neither bosonic nor
fermionic commutation relations, which leads to exciton-exciton Coulomb
interactions [32]. The excitations of the 2D material are, therefore,
described as interacting bosons with raising and lowering operators
$b^{\dagger}_{\mathbf{k}}$ and $b_{\mathbf{k}}$, respectively. We consider the
case of an electromagnetic resonator with a single resonance in the bandwidth
of interest for which the corresponding quantized electromagnetic field can be
written in terms of the photon raising and lowering operators $a^{\dagger}$
and $a$, respectively. In the framework of the ERC, the interaction between
excitons in the 2D material and the cavity mode is collected inside a single
exciton with annihilation operator
$B_{0}=(\sqrt{\sum_{\mathbf{k}}{|g_{\mathbf{k}}|^{2}}})^{-1}\sum_{\mathbf{k}}g^{*}_{\mathbf{k}}b_{\mathbf{k}}$,
where $g_{\mathbf{k}}$ is the coupling strength [19]. $B_{0}$ is, in turn,
coupled to the exciton reservoir described by the operators
$B_{i}=\sum_{\mathbf{k}}U_{i,\mathbf{k}}b_{\mathbf{k}}$, where
$U_{i,\mathbf{k}}$ is the transformation matrix. The first row of the ERC
transformation $U$, is
$U_{0,\mathbf{k}}=(\sqrt{\sum_{\mathbf{k}}{|g_{\mathbf{k}}|^{2}}})^{-1}\sum_{\mathbf{k}}g^{*}_{\mathbf{k}}$
while the remaining ones are constructed via Gram-Schmidt orthogonalisation as
orthonormal vectors such that $U$ is unitary. Furthermore, the part of the
Hamiltonian governing the modes with $i>0$ is re-diagonalized. The power of
the ERC thus relies on the fact that the exciton reservoir does not interact
directly with the cavity field. A formal introduction to the exciton reaction
coordinate in this kind of system is beyond the scope of this article, and we
refer to [19] for an in-depth discussion.

The Hamiltonian of the system describes the single isolated exciton and the
quantized electromagnetic field and can be written

$H^{\prime}_{\mathrm{s}}=\hbar\Omega_{0}B^{\dagger}_{0}B_{0}+\hbar\omega
a^{\dagger}a+\hbar G_{0}(B^{\dagger}_{0}a+a^{\dagger}B_{0}),$ (1)

where $\hbar$ is the reduced Planck constant, and the first two terms are the
free energy of the exciton reaction coordinate and the optical field,
respectively. $G_{0}$ is the light-matter coupling strength, which, in
practice, is obtained by integrating the optical mode distribution as detailed
in Ref. [19].

In this work, we are particularly interested in the so-called strong coupling
regime, where the coupling strength $G_{0}$ is so large that the exchange rate
between $B_{0}$ and $a$ is faster than the rate of any competing dissipative
process. In this regime, when the cavity energy is tuned across the exciton
energy, the dispersion relations for light and matter anticross, which
signifies the generation of two new states, known as the lower and the upper
polariton branches. These are quasi-particles made of both light and matter.
In the present model, the electromagnetic field is described by a single
optical mode inside the cavity, which implies that only a single upper and
lower polariton state is created. This is very different from the case of a
quantum well embedded between Bragg mirrors, where continuous polariton
branches arise [10].

To highlight the formation of the polariton quasiparticles, one can perform
the Hopfield-Bogoliubov transformation [20, 21] to recover the normal modes of
the Hamiltonian described by Eq. (1). These are

$\begin{split}&p_{\mathrm{lp}}=X_{\mathrm{lp}}B_{0}-C^{*}_{\mathrm{lp}}a\\\
&p_{\mathrm{up}}=C_{\mathrm{lp}}B_{0}+X^{*}_{\mathrm{lp}}a,\end{split}$ (2)

where $X_{\mathrm{lp}}$ and $C^{*}_{\mathrm{lp}}$ are the excitonic and
photonic weights of the new quasi-particles. Equation (1) can thus be written
in diagonal form

$H^{\prime}_{\mathrm{s}}=\epsilon_{\mathrm{lp}}p^{\dagger}_{\mathrm{lp}}p_{\mathrm{lp}}+\epsilon_{\mathrm{up}}p^{\dagger}_{\mathrm{up}}p_{\mathrm{up}},$
(3)

where $p^{\dagger}_{\mathrm{lp/up}}$ and $p_{\mathrm{lp/up}}$ are the raising
and lowering operators of the lower and upper polariton state, and
$\epsilon_{\mathrm{lp/up}}$ are the corresponding energies.

In order to account for the fact that excitons and polaritons are interacting
bosons, we add the leading interaction term describing Coulomb scattering to
the Hamiltonian in Eq. (1). This accounts for self-interaction within the
polariton state and reads [19]

$\hat{W}_{0}=W^{\prime}_{0}|X_{\mathrm{lp}}|^{4}p^{\dagger}_{\mathrm{lp}}p^{\dagger}_{\mathrm{lp}}p_{\mathrm{lp}}p_{\mathrm{lp}},$
(4)

with the effective non-linear interaction strength

$\footnotesize
W^{\prime}_{0}=\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}W_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}U_{0,\mathbf{k}+\mathbf{q}}U_{0,\mathbf{k}^{\prime}-\mathbf{q}}U^{*}_{0,\mathbf{k}^{\prime}}U^{*}_{0,\mathbf{k}}.$
(5)

From now on, the upper polariton branch will be neglected because it resides
high in energy and does not affect the dynamics described here. The lower
polariton will then be labelled by $p_{0}$ instead of $p_{\mathrm{lp}}$ to
highlight the relation with the ERC-operator $B_{0}$. See Eq. (2).

The Hamiltonian for the residual excitons $B_{i}$ forming the reservoir is

$H_{\mathrm{R}}=\hbar\sum_{i>0}\Omega_{i}B^{\dagger}_{i}B_{i}.$ (6)

The Coulomb interaction of the excitons is neglected here. This follows from
the assumption that the only region of the exciton Hilbert space, where the
exciton density is sufficiently large to give a significant contribution, is
the exciton-exciton interactions within the ERC, as described by Eq. (4). When
the characteristic confinement length of the electromagnetic field $L$, here
between $20\ \mathrm{nm}$ and $70\ \mathrm{nm}$, is large compared to the
exciton Bohr radius, $a_{\mathrm{B}}\approx 1.95\ \mathrm{nm}$, the
transformation elements $U_{0,\mathbf{k}}$ decay on a momentum scale that is
small compared to the momentum variation of
$W_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}$. In this regime it is
reasonable to use the approximation
$W_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}\simeq W_{000}$ [19, 33].

The polariton state interacts with the residual excitons of the environment
through Coulomb scattering. We consider here the primary interaction in which
two excitons in the reservoir scatter, generating a particle in the lower
polariton state and a higher-energy exciton. The opposite process, where a
polariton and a high-energy exciton generate two lower-energy excitons, is
also taken into account. Fig. 3 shows a schematic representation of such
processes between the continuum exciton reservoir and the single isolated
polariton state. The expression for $H_{\mathrm{int}}$ then reads

$H_{\mathrm{int}}=\sum_{j,i^{\prime},j^{\prime}}(\mathcal{U}_{j,i^{\prime},j^{\prime}}X_{\mathrm{lp}}p^{\dagger}_{0}B^{\dagger}_{j}B_{i^{\prime}}B_{j^{\prime}}+\mathrm{h.c.}),$
(7)

in which $\mathrm{h.c.}$ denotes the hermitian conjugate, and

$\footnotesize\mathcal{U}_{j,i^{\prime},j^{\prime}}=\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}W_{000}U_{0,\mathbf{k}+\mathbf{q}}U_{j,\mathbf{k}^{\prime}-\mathbf{q}}U^{*}_{i^{\prime},\mathbf{k}^{\prime}}U^{*}_{i,\mathbf{k}}.$
(8)

It follows that $W_{000}$ can be removed from the summation since it loses its
dependence on the exciton in-plane momentum. This approximation will be
beneficial later on.

Figure 4: Sketch of the scattering between two excitons in the reservoir.
After the scattering, a particle is captured into the lower polariton state,
and an exciton is scattered to a high-energy state in the reservoir. The
subscript $i$ labels the bosonic states in the ERC and replaces the in-plane
momentum $\mathbf{k}$ in the new basis.

###  B Born-Markov Master Equation

In the formalism of open quantum systems, the 2D material coupled to the
nanocavity can be modelled with the Hamiltonian

$H=H_{\mathrm{s}}+H_{\mathrm{R}}+H_{\mathrm{int}}$ (9)

in which the first term describes the system ERC-mode Hamiltonian in the
polariton basis, in which $H_{\mathrm{s}}=H^{\prime}_{\mathrm{s}}+W_{0}$ as
defined in Eqs. (3) and (4); $H_{\mathrm{R}}$ is given by Eq. (6) and encodes
the exciton reservoir not interacting with the cavity mode; lastly,
$H_{\mathrm{int}}$ is the system-environment interaction term described by Eq.
(7).

The second-order perturbative Born-Markov master equation in the Schröedinger
picture, obtained by tracing out the environment, is given by [34]

$\footnotesize\frac{\partial\rho_{\mathrm{s}}(t)}{\partial
t}=-\frac{i}{\hbar}[H_{\mathrm{s}},\rho_{\mathrm{s}}(t)]-\frac{1}{\hbar^{2}}\int^{\infty}_{0}\mathrm{d}\tau\mathrm{Tr_{\mathrm{R}}}[H_{\mathrm{int}},[\tilde{H}_{\mathrm{int}}(-\tau),\chi(t)]],$
(10)

where $\rho_{\mathrm{s}}$ is the system density matrix. The integrand is
traced over the reservoir, and $\tilde{H}_{\mathrm{int}}(-\tau)$ is a time-
dependent Hamiltonian operator. Eq. (10) has the explicit expression

$\footnotesize\begin{split}\frac{\partial\rho_{\mathrm{s}}(t)}{\partial
t}=&-\frac{\mathrm{i}}{\hbar}[H_{\mathrm{s}},\rho_{\mathrm{s}}(t)]+\frac{1}{\hbar^{2}}\int^{\infty}_{0}\mathrm{d}\tau\Big{(}(p^{\prime}_{0}\rho_{\mathrm{s}}(t)p^{\dagger}_{0}-p^{\dagger}_{0}p^{\prime}_{0}\rho_{\mathrm{s}}(t))\left<F^{\dagger}_{0}F^{\prime}_{0}\right>\\\
&+(p^{\prime\dagger}_{0}\rho_{\mathrm{s}}(t)p_{0}-p_{0}p^{\prime\dagger}_{0}\rho_{\mathrm{s}}(t))\left<F_{0}F_{0}^{{}^{\prime}\dagger}\right>-\mathrm{h.c.}\Big{)},\end{split}$
(11)

where the operators labelled with prime are still in the interaction picture
while

$\footnotesize
F^{\dagger}_{0}=\sum_{ijh}\mathcal{U}_{i,j,h}B^{\dagger}_{i}B_{j}B_{h}$ (12)

describes the scattered excitons, Eq. (8).

###  C Rate equations

To study the dynamics of the lower polariton population
$N_{0}=\left<p^{\dagger}_{0}p_{0}\right>$, we follow Ref. [35] to develop a
set of rate equations. We generalize this model to take into account the
breaking of the in-plane translational symmetry induced by the nanocavity
[19]. The evolution of the exciton density $n_{x}$ and the reservoir
temperature $T_{x}$ are also described, leading to a set of self-consistent
rate equations.

The time evolution of the expectation value of the lower polariton population
is given as
$\dfrac{\mathrm{d}}{\mathrm{d}t}\left<p^{\dagger}_{0}p_{0}\right>=\mathrm{Tr}\Big{(}p^{\dagger}_{0}p_{0}\dfrac{\partial\rho_{\mathrm{s}}}{\partial_{t}}\Big{)}$,
which leads to the equation

$\footnotesize\begin{split}\frac{dN_{0}}{\mathrm{d}t}=&-\frac{1}{\hbar^{2}}\int^{\infty}_{0}\mathrm{d}\tau\Big{(}\left<p^{\dagger}_{0}p_{0}p^{\dagger}_{0}p^{\prime}_{0}\right>-\left<p^{\dagger}_{0}p^{\dagger}_{0}p_{0}p^{\prime}_{0}\right>\Big{)}\left<F^{\dagger}_{0}F^{\prime}_{0}\right>\\\
&+\Big{(}\left<p^{\dagger}_{0}p_{0}p_{0}p^{{}^{\prime}\dagger}_{0}\right>-\left<p_{0}p^{\dagger}_{0}p_{0}p^{{}^{\prime}\dagger}_{0}\right>\Big{)}\left<F_{0}F^{{}^{\prime}\dagger}_{0}\right>+\mathrm{h.c.}.\end{split}$
(13)

The multiple operator expectation values are expanded using Wick’s theorem
[36]. For the six operators in $\left<F^{\dagger}_{0}F^{\prime}_{0}\right>$
and $\left<F_{0}F^{{}^{\prime}\dagger}_{0}\right>$, the assumption of a
thermalized Maxwell-Boltzmann reservoir allows us to simplify the computation.
We make the assumption that the electromagnetic mode function is both
separable and localized. In this approximation, the precise overlap of the
exciton reaction coordinate and the resonant field results in a coupling
strength $G_{0}$ that remains unaffected by the electromagnetic field’s
lateral extent [19]. Furthermore, the lateral field distribution of the
electromagnetic field is taken as Gaussian, which simplifies the computation
of the Coulomb in-scattering rate and the out-scattering rate of the lower
polariton state. The full derivation is reported in Appendix A. The final
expression for the time-dependence of $N_{0}$ then reads

$\begin{split}\frac{dN_{0}}{dt}=&(N_{0}+1)n^{2}_{x}W_{\mathrm{in}}-N_{0}n_{x}W_{\mathrm{out}_{\mathrm{1}}}-N_{0}n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}\\\
&+[n^{3}_{x}W_{\mathrm{in_{c}}}-n^{2}_{x}N_{0}W_{\mathrm{out_{c}}}]-\Gamma_{\mathrm{lp}}N_{0},\end{split}$
(14)

where $W_{\mathrm{in}}$ and $W_{\mathrm{in}_{\mathrm{c}}}$ are the two rates
of scattering into the lower polariton, scaling as $n_{x}^{2}$ and
$n_{x}^{3}$, respectively. $W_{\mathrm{out}_{\mathrm{1}}}$,
$W_{\mathrm{out}_{\mathrm{2}}}$ and $W_{\mathrm{out}_{\mathrm{c}}}$ are the
out-scattering rates, the former rate scaling with $n_{x}$ and the others with
$n_{x}^{2}$. In the following, we use the label $W$ to refer to all the
scattering rates. Crucially, all the $W$-coefficients depend on the exciton
temperature. The rates are thus not fixed quantities but are time-dependent.
The last term of Eq. (14) has been added phenomenologically to account for the
radiative losses of the polaritons. In particular,
$\Gamma_{\mathrm{lp}}=C^{2}_{\mathrm{lp}}/\tau_{\mathrm{ph}}+X^{2}_{\mathrm{lp}}/\tau_{x}$
where $\tau_{\mathrm{ph}}$ is the cavity photon lifetime and
$\tau_{x}=\Gamma^{-1}_{x}$ is the exciton lifetime [37]. Using the Maxwell-
Boltzmann distribution, we define the exciton density $n_{x}$ through

$N^{x}=\frac{2\pi n_{x}}{\rho_{x}k_{\mathrm{B}}T_{x}},$ (15)

in which $N^{x}$ is the total occupation number 111the $i_{th}$ exciton state
occupation number is $N^{x}_{i}=N^{x}\exp{-\beta\epsilon_{i}}$, with
$\beta=(T_{x}k_{\mathrm{B}})^{-1}$ and $\epsilon_{i}$ the corresponding
energy., $\rho_{x}=(m_{\mathrm{h}}+m_{\mathrm{e}})/\hbar^{2}$ is the 2D
density of state that does not depend on space but only on the electron and
hole masses $m_{\mathrm{e}}$ and $m_{\mathrm{h}}$, respectively. Finally,
$k_{\mathrm{B}}$ is the Boltzmann constant, and $T_{x}$ is the exciton
reservoir temperature, which evolves with the pumping as discussed below.

A self-consistent model for the dynamics requires the inclusion of two
additional equations. These govern the changes in the density and energy of
the exciton reservoir and are derived from Eq. (14) by accounting for
particles added to or removed from the system. The evolution of the exciton
density involves a non-resonant pump to account for the injection of
thermalized excitons into the reservoir, $p_{x}$. Whenever a particle
undergoes scattering into the lower polariton state, an exciton is extracted
from the reservoir, and a second exciton with higher energy is produced, as
depicted in Fig. 4. This process results in increasing the reservoir
temperature. The third rate accounts for these dynamical changes in the
exciton reservoir temperature. Consequently, the dynamics of $n_{x}$ can be
described as follows

$\begin{split}\footnotesize\frac{dn_{x}}{dt}&=-(N_{0}+1)n^{2}_{x}\frac{W_{\mathrm{in}}}{S}+N_{0}n_{x}\frac{W_{\mathrm{out}_{\mathrm{1}}}}{S}+\\\
&N_{0}n^{2}_{x}\frac{W_{\mathrm{out}_{\mathrm{2}}}}{S}+p_{x}-n_{x}\Gamma_{x}+N_{0}n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}-\\\
&[n^{3}_{x}W_{\mathrm{in_{c}}}-n^{2}_{x}N_{0}W_{\mathrm{out_{c}}}]\frac{1}{S}\end{split}$
(16)

where $S$ is the area of the 2D sheet. The non-resonant pump rate $p_{x}$ is
modelled as a flux, having dimensions of inverse time per inverse square meter
and scaling with the area of the 2D sheet. In the absence of polaritons, the
density of the excitons in the reservoir created by the pump is then
$n_{x}=p_{x}\tau_{x}$. For $T_{x}=e_{x}/k_{\mathrm{B}}$

$\begin{split}\footnotesize\frac{de_{x}}{dt}&=-(N_{0}+1)n^{2}_{x}\frac{W_{\mathrm{in}}}{S}\epsilon_{\mathrm{lp}}+N_{0}n_{x}\frac{W_{\mathrm{out}_{\mathrm{1}}}}{S}\epsilon_{\mathrm{lp}}+\\\
&N_{0}n^{2}_{x}\frac{W_{\mathrm{out}_{\mathrm{2}}}}{S}\epsilon_{\mathrm{lp}}+p_{x}k_{\mathrm{B}}T_{\mathrm{L}}-n_{x}\Gamma_{x}k_{\mathrm{B}}T_{x}-\\\
&[n^{3}_{x}W_{\mathrm{in_{c}}}-n^{2}_{x}N_{0}W_{\mathrm{out_{c}}}]\frac{\epsilon_{\mathrm{lp}}}{S}\end{split}$
(17)

where $e_{x}$ is the energy density of the exciton reservoir, which evolves
according to the equation, while $\epsilon_{\mathrm{lp}}$ is the energy of the
single-state lower polariton branch. The term
$p_{x}k_{\mathrm{B}}T_{\mathrm{L}}$, is considered here to be given by the
injection of excitons at the lattice temperature $T_{\mathrm{L}}$. One should
stress the difference between the lattice temperature, which accounts for the
temperature at which the experiment is performed, and the exciton temperature,
which is a dynamical variable.

Aside from the inclusion of the second-order out-scattering rate,
$W_{\mathrm{out}_{\mathrm{2}}}$, Eqs. (14), (16) and (17) have two extra terms
compared to previous models developed for quantum wells using the same methods
[35]. Specifically, two additional scattering rates, $W_{\mathrm{in_{c}}}$ and
$W_{\mathrm{out_{c}}}$, arise due to the loss of translation invariance in the
system, see Appendix A.

###  D Scattering rates

Here, we discuss the physics of the Coulomb scattering rates entering Eqs.
(14), (16) and (17). We also explain the impact of the loss of translation
invariance on the system’s dynamics.

Due to the translational symmetry breaking the scattering rates do not have an
analytic expression but must be computed numerically. We report here the
expression for $W_{\mathrm{in}}$. Further details can be found in Appendix B
together with expressions for the other rates,

$\begin{split}W_{\mathrm{in}}=&\frac{4\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}\rho_{x}}{2\pi(k_{\mathrm{B}}T_{x})^{2}}\\\
&\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}\mathrm{d}\theta_{2}\mathrm{d}\theta_{3}e^{-2L^{2}(k^{\prime}_{1}k_{2}cos\theta_{2}-k^{\prime}_{1}k_{3}cos\theta_{3}}e^{-L^{2}k^{{}^{\prime}2}_{1}}\\\
&e^{-\beta(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}})}\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}).\end{split}$
(18)

The integrals are performed over the exciton energies
$\epsilon_{\mathbf{k}_{2},\mathbf{k}_{3}}$ and the angles $\theta_{2,3}$;
$\mathcal{H}$ is the Heaviside step function, and $k^{\prime
2}=k^{2}_{3}-k^{2}_{2}+\frac{2m\epsilon_{\mathrm{lp}}}{\hbar^{2}}$. To reduce
the number of integrals we use the approximation introduced in Ref. [19],
which assumes the electromagnetic field to be separable with a Gaussian
lateral field distribution. This simplifies the expression for the light-
matter coupling $G_{0}$, and enables to find an analytic expression for the
ERC transformation matrix
$|U_{0,\mathbf{k}}|^{2}\approx\frac{L^{2}}{S}e^{-L^{2}k^{2}}$, see Appendix B.

Figure 5: Numerically computed in-scattering rate for $G_{0}=3.8\mathrm{meV}$
and different optical cavity in-plane confinements. The black dotted line is
the result expected for a quantum-well model where translation invariance is
preserved, in which case it is possible to derive an analytic expression for
$W_{\mathrm{in}}$ [10]. The inset shows the behaviour of the scattering-in
rate for $L=20\ \mathrm{nm}$, for three different coupling strengths.

The example in Fig. 2 shows a so-called bowtie cavity displaying extreme
dielectric confinement for which the field distribution of the quasi-normal
mode can be found in Ref. [12]. We note here that the description of the
transverse spatial profile of the field distribution by a Gaussian function is
not necessarily a good approximation. Nevertheless, the analysis suggests that
$G_{0}$ depends chiefly on the characteristic out-of-plane confinement length
of the cavity, $L^{-1}_{z}$, whereas
$W_{0}=W^{\prime}_{0}|X_{\mathrm{lp}}|^{4}$ depends mostly on the in-plane
confinement length, $L$, and on none of the other length scales of the system
[19]. In principle, this enables tuning the light-matter coupling strength
just by changing the value of $L_{z}$. The three lengths scale $L$, $L_{z}$,
$S$ that enter $W$, then play specific roles that are well distinguished.

Eq. (18) shows the explicit dependence of the rates on the reservoir
temperature. This is plotted in Fig. 5 for different in-plane confinement
lengths of the optical cavity. The scattering rates depend on the exciton
occupation number $N^{x}_{i}$, which follows Maxwell-Boltzmann statistics. The
dynamics of Eq. (17) then also determines the magnitude of the rates. The
black dotted line in Fig. 5 shows the analytic expression for
$W_{\mathrm{in}}$ found in Ref. [10], where the spatial confinement of the
electromagnetic field is not taken into account, and it is assumed that the
surface confinement is equal to $S$. Such analytic expression has the same
temperature dependence obtained from our numerical simulation. Finally, the
observed increase in scattering rates with varying in-plane cavity lengths can
be attributed to the broadening of the Gaussian mode profile. This broadening
results in an increased overlap with the Maxwell-Boltzmann distribution of the
reservoir excitons.

Figure 6: Sketch of the scattering between a lower polariton and a higher
energy exciton. After the scattering, two lower exciton states are filled. The
lower polaritons can scatter with any exciton in the reservoir as long as
energy is conserved in the process.

As previously mentioned, the lack of conservation of in-plane momentum implies
that in the expansion for $\left<F^{\dagger}_{0}F_{0}\right>$ and
$\left<F_{0}F^{\dagger}_{0}\right>$, extra terms appear that would normally
cancel out for translationally invariant systems. In the following, these are
labelled by $W_{\mathrm{in}_{\mathrm{c}}}$ and
$W_{\mathrm{out}_{\mathrm{c}}}$; see Appendix A and B. Physically, we can
understand these terms as an increase in the value of the scattering rates due
to the presence of a single isolated state rather than a continuous polariton
branch. In a translationally invariant 2D system, such as a quantum well in a
microcavity, Coulomb interactions conserve both energy and momentum. Two
excitons in the reservoir, denoted by momentum and energy $k^{x}_{\mathrm{i}}$
and $\epsilon^{x}_{i}$ respectively, scatter into a lower polariton state
$(k^{\mathrm{lp}}_{\mathrm{f}},\epsilon^{\mathrm{lp}}_{k})$ and a higher-
energy exciton $(k^{x}_{\mathrm{f}},\epsilon^{x}_{\mathrm{f}})$. Considering
that the wave vector of the lower polariton is negligible compared to the
scattering excitons in the reservoir, $k^{\mathrm{lp}}_{\mathrm{f}}\approx 0$,
the conservation of momentum and energy results in the following conditions:

$\begin{split}&2|k^{x}_{\mathrm{i}}|=|k^{x}_{\mathrm{f}}|\rightarrow
4\epsilon^{x}_{\mathrm{i}}=\epsilon^{x}_{\mathrm{f}}\\\
&2\epsilon^{x}_{\mathrm{i}}=-|\epsilon^{\mathrm{lp}}_{k}|+\epsilon^{x}_{\mathrm{f}}\rightarrow\epsilon^{x}_{\mathrm{i}}=|\epsilon^{\mathrm{lp}}_{k}|/2,\epsilon^{x}_{\mathrm{f}}=2|\epsilon^{\mathrm{lp}}_{k}|\end{split}$
(19)

This implies that a minimum
$\epsilon^{x}_{\mathrm{i}}=\epsilon^{x}_{\mathrm{lp}}/2$ is required for the
scattering processes described, where now $\epsilon^{x}_{\mathrm{lp}}$ labels
one of the states of the lower polariton branch, see Ref. [10]. These
considerations do not apply to the system we study, where there is no
condition on the initial energy of the scattering excitons
$\epsilon^{x}_{\mathrm{i}}$. As illustrated in Fig. 6, a polariton can
interact with any exciton in the reservoir, and any exciton can scatter into
the lower isolated state, provided that energy is conserved.

’

Figure 7: Numerically computed $W_{\mathrm{in}_{\mathrm{c}}}$ scattering rate
for $G_{0}=3.8\ \mathrm{meV}$ and different optical cavity in-plane
confinement lengths. Increasing $L$, $W_{\mathrm{in}_{\mathrm{c}}}$ initially
grows faster, but increasing the in-plane cavity length even more, its value
remains stable.

Fig. 7 shows how the correction to the in-scattering rate varies with
temperature for different in-plane confinement values. As the cavity size
increases, the rate follows the same trend as $W_{\mathrm{in}}$. Furthermore,
its value remains stable for $L=50\ \mathrm{nm},70\ \mathrm{nm}$. It follows
that when $L$ becomes sufficiently large, the correction terms cease to grow
and become less relevant in the system’s dynamics.

The validity of the model is subject to two important constraints. The first
one pertains to the maximum allowable value of the exciton density, which is
limited to less than the saturation density $n_{\mathrm{sat}}\approx
1/a^{2}_{\mathrm{B}}=2\times 10^{17}\ \mathrm{m}^{-2}$. This restriction
ensures that the system remains below the Mott transition and justifies the
description of the excitons as interacting bosons [9, 39]. The second
constraint limits the maximum allowable value of the total exciton occupation
number. The population of the exciton reservoir is modelled as a Maxwell-
Boltzmann distribution where quantum degeneracy effects are negligible. This
requires $N^{x}<1$ and from Eq. (15) implies that
$n_{x}<k_{\mathrm{B}}T_{x}\rho_{x}/2\pi$ [10].

##  III Results

In this section, we present and discuss numerical results obtained from our
model. In the simulations below, the 2D material is chosen to be Tungsten
disulfide ($\mathrm{WS_{2}}$), for which the important material parameters are
shown in Table 1. For the nanocavity, the out-of-plane-confinement of the
model-cavity has been adjusted to tune the value of the coupling strength
$G_{0}$ [19].

$a_{\mathrm{B}}[\mathrm{nm}]$ | $m_{\mathrm{h}}[m_{0}]$ | $m_{\mathrm{e}}[m_{0}]$ | $E_{\mathrm{g}}[\mathrm{eV}]$ | $E_{\mathrm{b}}[\mathrm{eV}]$ | $\tau_{x}[\mathrm{ps}]$
---|---|---|---|---|---
$1.95$ | $0.34$ | $0.33$ | $2.53$ | $0.52$ | $2$
Table 1: Relevant material parameters for a monolayer, $\mathrm{WS_{2}}$. The
listed quantities are the exciton Bohr radius $a_{\mathrm{B}}$; the effective
hole and electron masses $m_{\mathrm{h}}$ and $m_{\mathrm{e}}$, respectively,
given in units of the free electron mass $m_{0}$; the band-gap energy
$E_{\mathrm{g}}$, and the exciton binding energy $E_{\mathrm{b}}$ [28, 40].
For the exciton lifetime $\tau_{x}$ we refer to Ref. [41].

We study and discuss the system at steady-state as a function of the pumping
and find that the evolution of the ground-state occupation exhibits a distinct
kink at which it can increase by orders of magnitude with only a moderate
increase in the pump rate. We regard this highly non-linear behaviour as the
onset of polariton condensation. We pay particular attention to the role
played by the temperature, differentiating the role of $T_{x}$ and
$T_{\mathrm{L}}$. The dynamics of the former turn out to be strongly linked to
the number of particles in the lower polariton state. By tuning the latter, we
can control the onset of condensation as we vary the pumping rate.

In steady-state, from Eqs. (14), (16) and (17), we get the following coupled
equations

$\displaystyle
N_{0}^{\mathrm{ss}}=\frac{R_{\mathrm{sp}}(n_{x}^{\mathrm{ss}},T_{x}^{\mathrm{ss}})}{\Gamma_{\mathrm{lp}}-G_{\mathrm{st}}(n_{x}^{\mathrm{ss}},T_{x}^{\mathrm{ss}})},$
(20) $\displaystyle
p_{x}=\Gamma_{x}n_{x}^{\mathrm{ss}}+\frac{\Gamma_{\mathrm{lp}}N_{0}^{\mathrm{ss}}}{S},$
(21) $\displaystyle
p_{x}k_{B}T_{\mathrm{L}}=\Gamma_{x}n_{x}^{\mathrm{ss}}k_{B}T_{x}^{\mathrm{ss}}+\frac{\Gamma_{\mathrm{lp}}N_{0}^{\mathrm{ss}}\epsilon_{\mathrm{lp}}}{S},$
(22)

where $R_{\mathrm{sp}}(n_{x},T_{x})$ and $G_{\mathrm{st}}(n_{x},T_{x})$ are
functions representing respectively the rate of spontaneous and stimulated
scattering into the polariton branch,

$\displaystyle R_{\mathrm{sp}}$
$\displaystyle=n_{x}^{2}W_{\mathrm{in}}(T_{x})+n_{x}^{3}W_{\mathrm{in}_{\mathrm{c}}}(T_{x}),$
(23)
$\displaystyle\begin{split}G_{\mathrm{st}}&=n_{x}^{2}W_{\mathrm{in}}(T_{x})-n_{x}^{2}W_{\mathrm{out}_{\mathrm{2}}}(T_{x})-\\\
&-n_{x}^{2}W_{\mathrm{out}_{\mathrm{c}}}(T_{x})-n_{x}W_{\mathrm{out}_{\mathrm{1}}}(T_{x}).\end{split}$
(24)

Eq. (20) gives the polariton number in steady-state for a given temperature
and exciton density. Eq. (21) represents the conservation of particles, and
Eq. (22) corresponds to the conservation of energy. The input-output relations
are found by solving these coupled nonlinear algebraic equations numerically.

Figure 8: Upper plot, polariton occupation as a function of pump power for
different in-plane confinement lengths of the optical cavity, $G_{0}=3.8\
\mathrm{meV}$ and $T_{\mathrm{L}}=8\ \mathrm{K}$. The small dots in each curve
label the point where the corresponding occupation number of the excitons
$N^{x}\sim 1$ and Maxwell-Boltzmann distribution breaks down in this limit.
The black star refers to the value corresponding to the Mott transition of
$n_{x}$. The black dotted line represents the system’s dynamics without
correction terms, assuming analytic expressions for in-scattering and out-
scattering rates [10] and considering a single polariton state. The bottom
plot shows the variation in the exciton reservoir temperature.

Fig. 8 shows the change in the steady-state values of the polariton occupation
(upper plot) and the reservoir temperature (bottom plot) as a function of the
input power. The four curves refer to a system with light-matter coupling
$G_{0}=3.8\ \mathrm{meV}$, $T_{\mathrm{L}}=8\ \mathrm{K}$ and different in-
plane lengths $L$ of the optical cavity. The extension of the lattice surface
is maintained constant at $S=(15\ \mathrm{\mu m})^{2}$ for the four
simulations so that the same pump is required to excite the same initial
density of excitons. $L$ is progressively increased, enhancing the value of
the rates as shown in Fig. 5 for $W_{\mathrm{in}}$ and Fig. 7 for
$W_{\mathrm{in_{c}}}$. The upper plot in Fig. 8 shows that as $L$ increases, a
higher pumping rate becomes necessary to observe a kink in the evolution of
$N_{0}$. This continues until $L=70\ \mathrm{nm}$ (red curve) where the
required pump decreases, bringing the curve closer to the black dotted one,
which represents the system’s dynamics for a single polariton state, without
correction terms and assuming analytic expressions for the in-scattering rate
(black dotted curve in Fig. 5) and the out-scattering rate [10]. We can
conclude that larger cavity sizes necessitate bigger pumping to observe a kink
in $N_{0}$. For $L$ sufficiently large, moreover, the effect of
$W_{\mathrm{in_{c}}}$ and $W_{\mathrm{out_{c}}}$ becomes progressively smaller
in $N_{0}$ dynamics and we recover the limit where $W_{\mathrm{in}}$ and
$W_{\mathrm{out_{1}}}$ have analytic expressions.

As seen in the bottom panel of Fig. 8, the onset of condensation leads to an
increase in the reservoir temperature because more excitons are scattered into
higher energy states. This increase in the value of $T_{x}$ can, in turn,
result in bistability in the evolution of $N_{0}$, as evident from the
dynamics displayed by the red and black dotted curves in the top panel Fig. 8.
Further elucidation of this phenomenon will be presented in subsequent
discussions.

Figure 9: Upper plot, the evolution of the scattering rates as a function of
the pump power $p_{x}$ for $G_{0}=3.8\ \mathrm{meV}$ and $L=20\ \mathrm{nm}$.
Middle plot, evolution of the total scattering-out rate
$\Gamma_{\mathrm{out}}(p_{x})=n^{2}_{x}W_{\mathrm{out}_{\mathrm{1}}}+n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}+n^{2}_{x}W_{\mathrm{out_{c}}}$
compared to
$n_{x}W_{\mathrm{out}_{\mathrm{1}}}+n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}$.
Bottom plot, corresponding evolution of the exciton density $n_{x}$. The
dotted black line is the $p_{x}$ value corresponding to the $N_{0}$ threshold.
The black stars indicate the onset of the Mott transition.

The evolution of $T_{x}$ illustrates that the system behaves qualitatively
differently from the traditional case of a system with translational
invariance. In the latter case, the exciton temperature increases with the
pumping [10]. From Fig. 8, however, the exciton temperature approaches a
saturation limit, which depends on the coupling strength and the cavity size.
This behaviour seemed rather linked to the consideration in Sec. D, where it
was pointed out that there is not a minimum energy required for the particle
to scatter into $p_{0}$.

Figure 10: Polariton occupation number $N_{0}$ as a function of the pump power
$p_{x}$ for different coupling strengths $G_{0}$ and lattice temperatures. In
the upper plot $T_{\mathrm{L}}=1\ \mathrm{K}$, in the middle
$T_{\mathrm{L}}=8\ \mathrm{K}$ and in the bottom $T_{\mathrm{L}}=12\
\mathrm{K}$. The small dots in each curve label the point where the
corresponding occupation number of the excitons $N_{x}\sim 1$ and Maxwell-
Boltzmann distribution breaks down. The black star indicates the density at
which the Mott transition for the exciton reservoir occurs. Dashed regions of
the curves indicate an exciton density higher than the limit imposed by
Maxwell-Boltzmann when the exciton energy is still defined by the initial
lattice temperature $T_{\mathrm{L}}$.

The different scattering rates, together with the evolution of $n_{x}$ as a
function of the pumping rate, are plotted in Fig. 9, for a cavity length of
$L=20\ \mathrm{nm}$. The black vertical line shows the pumping rate needed to
reach the threshold. Here, we define the threshold as the condition for which,
for a given temperature, the gain equals the losses, see Eq. (29) below. From
the top plot of Fig. 9, we observe that $W_{\mathrm{in_{c}}}n^{3}_{x}$ grows
progressively faster than the other rates as a function of the pump. The total
Coulomb out-scattering of the polariton state is also modified by
$W_{\mathrm{out}_{\mathrm{c}}}$, as seen in the centre plot. In Sec. II, we
interpret the emergence of the correction terms as an enhancement of the
overall scattering rates originating from a single polariton state; see Fig.
6. To further illustrate this behaviour, Fig. 9 shows the evolution of
$\Gamma_{\mathrm{out}}(p_{x})=n_{x}W_{\mathrm{out}_{\mathrm{1}}}+n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}+n^{2}_{x}W_{\mathrm{out_{c}}}$
compared to
$n_{x}W_{\mathrm{out}_{\mathrm{1}}}+n^{2}_{x}W_{\mathrm{out}_{\mathrm{2}}}$.
From these comparisons, it becomes evident that the correction rates, in this
case $n_{x}^{3}W_{\mathrm{out_{c}}}$, become increasingly more important as
the system approaches the threshold, emphasizing the necessity of accounting
for the loss of translation invariance to describe the system accurately. We
therefore conclude that we cannot neglect these terms since they contribute in
a non-negligible way to the dynamics.

To gain deeper insights into the dynamics, we compare the steady-state values
of the lower polariton population $N_{0}$ for three different coupling
strengths while keeping the confinement length fixed at $20\ \mathrm{nm}$, see
Fig. 10. This analysis emphasizes the crucial role played by the lattice
temperature in defining a kink and the presence of hysteresis. In the upper
plot of the figure, in which we consider the case in which the initial
temperature is set to $T_{\mathrm{L}}=1\ \mathrm{K}$, a distinct kink is
observed. As we progressively increase $T_{\mathrm{L}}$, however, the middle
and bottom plots show that only the systems with higher values of $G_{0}$
exhibit a deviation from linear evolution. In the present model, as well as in
a translation-invariant case, a higher coupling strength, meaning an increase
in the value of $|\epsilon_{\mathrm{lp}}|$, leads to the scattering of
reservoir excitons to higher-energy states, raising the reservoir temperature
even more. Based on this understanding, we can make some conclusions on the
saturation temperature earlier introduced in Fig. 8. This value is expected to
grow with $G_{0}$ as well. Thus, we expect to observe the same dynamics at
higher temperatures as $G_{0}$ increases.

Figure 11: Colour map of the steady-state value of the ground state
population, $N_{0}^{\mathrm{ss}}$ versus the exciton reservoir density $n_{x}$
and temperature $T_{x}$. The black line shows the threshold exciton density
$n_{x,\mathrm{th}}(T_{x})$, where $N_{0}^{\mathrm{ss}}$ diverges. The red
line, $T_{\mathrm{l}}\simeq 12\ \mathrm{K}$, defines the limit temperature,
behind which the system cannot reach a threshold anymore.

A formal definition for the limit temperature $T_{\mathrm{l}}$ can be obtained
by studying the polariton distribution at steady-state. Equation (14) can be
rewritten as

$\frac{dN_{0}}{dt}=R-\Gamma_{\mathrm{lp}}N_{0},$ (25)

where the polariton generation rate $R$ containes a spontaneous contribution
and a stimulated contribution

$R=R_{\mathrm{sp}}+G_{\mathrm{st}}N_{0}.$ (26)

both depending on $T_{x}$ and $n_{x}$ and previously defined in Eqs. (23) and
(24). In steady-state, if $n_{x}$ and $T_{x}$ are known, the polariton
population is given by

$N_{0}^{\mathrm{ss}}=\frac{R_{\mathrm{sp}}}{\Gamma_{\mathrm{lp}}-G_{\mathrm{st}}},$
(27)

which is plotted in Fig. 11 as a function of $n_{x}$ and $T_{x}$. For small
pumping, we can approximate

$N_{0}^{\mathrm{ss}}\approx R_{\mathrm{sp}}/\Gamma_{\mathrm{lp}}.$ (28)

This approximation breaks down near the threshold. For a given temperature
$T_{x}$, we define the threshold exciton density $n_{x,\mathrm{th}}$ as the
point where the gain equals the losses

$G_{\mathrm{st}}(n_{x,\mathrm{th}}(T_{x}),T_{x})-\Gamma_{\mathrm{lp}}=0,$ (29)

which is plotted as the black line in Fig. 11. Solving Eq. (29) for
$n_{x,\mathrm{th}}$, we get the solution

$\footnotesize
n_{x,\mathrm{th}}=\frac{W_{\mathrm{out}_{\mathrm{1}}}}{2W_{2}}+\sqrt{\left(\frac{W_{\mathrm{out}_{\mathrm{1}}}}{2W_{2}}\right)^{2}+\frac{\Gamma_{\mathrm{lp}}}{W_{2}}},$
(30)

with the requirement for a threshold, i.e. a positive and real solution, being
$W_{2}=W_{\mathrm{in}}-W_{\mathrm{out}_{\mathrm{2}}}-W_{\mathrm{out}_{\mathrm{c}}}>0$.
At $T_{x}=T_{\mathrm{l}}$ the threshold diverges, coinciding with the point
where $W_{2}\rightarrow 0$, which is shown as the red line in Fig. 11. Above
this temperature, $G_{\mathrm{st}}$ is negative for all $n_{x}$ and we do not
obtain stimulated emission into the polariton branch. Thus, $T_{\mathrm{l}}$
is defined by

$W_{2}(T_{\mathrm{l}})\equiv 0.$ (31)

Figure 12: Three-dimensional plot of the steady-state values of $N_{0}$ versus
$p_{x}$ and $T_{\mathrm{L}}$, for fixed $G_{0}=3.8\ \mathrm{meV}$ and $L=20\
\mathrm{nm}$. The colour map represents the variation in the polariton
occupation number, lighter colours correspond to higher values of $N_{0}$.
There are three phases corresponding to bistability, threshold and not-
threshold. These are reached by increasing the value of the lattice
temperature $T_{\mathrm{L}}$ and letting $N_{0}$ evolve as a function of the
pumping $p_{x}$.

This result can be recognized in the three-dimensional plot in Fig. 12, which
shows the steady-state values of $N_{0}$ versus $p_{x}$ and $T_{\mathrm{L}}$,
for fixed $G_{0}=3.8\ \mathrm{meV}$ and $L=20\ \mathrm{nm}$. For smaller
values of $T_{\mathrm{L}}$, the dynamics of $N_{0}$ show initial bistability
followed only by a kink. Increasing the temperature above
$T_{\mathrm{L}}\simeq 12\ \mathrm{K}$, even the kink in $N_{0}$ disappears,
and its growth as a function of the pump becomes linear in a log-log scale. We
identify $T_{\mathrm{l}}$ as the value of the temperature marking the boundary
between these two regions.

###  A Discussion

The simulations show that a kink in the polariton population is not always
observed, even when the occupation number reaches values larger than one, and
we expect final state stimulation as a consequence of the bosonic nature of
the quasi-particles [42]. The stimulated scattering into the ground state
depends on the initial lattice temperature, and above some limit value
$T_{\mathrm{l}}$, the dynamics do not show any threshold, as is visible in
Fig. 12. We also point out that for a sufficiently low initial temperature,
the exciton reservoir evolves to $T_{\mathrm{l}}$ and maintains this
temperature regardless of how much the pumping rate is increased. As noted
above, this is a consequence of the loss of translation invariance in the
system, leading to a result different from the one predicted for quantum wells
[10].

When the initial lattice temperature is further decreased, Fig. 12 shows the
appearance of hysteresis. Bistability is often observed in polariton systems
that are resonantly pumped [43]. This originates from a nonlinear frequency
shift of the polariton mode along with near-resonant, but slightly detuned,
pumping. Initially, the threshold is high because of the detuning, but if a
polariton population is built up, the polariton energy shifts into resonance
with the pump, decreasing the threshold [44]. We consider a non-resonant pump.
Also, the magnitude of $G_{0}$, which defines $|\epsilon_{\mathrm{lp}}|$, is
fixed in the beginning by $L_{z}$ [19], so the energy of the lower polariton
state is not shifting. Instead, the presence of bistability appears to be
related to the dynamics of $T_{x}$. In particular, the dependence of the
scattering rates on the exciton temperature induces a temperature-dependent
threshold $n_{x,\mathrm{th}}(T_{x})$. Bistability occurs if
$n_{x,\mathrm{th}}(T_{x})$ decreases as $T_{x}$ increases above the initial
lattice temperature. This becomes more clear from Fig. 11 where bistability
can be observed if $T_{\mathrm{L}}$ is chosen such that there exists a region
of $T_{x}>T_{\mathrm{L}}$ where
$n_{x,\mathrm{th}}(T_{x})<n_{x,\mathrm{th}}(T_{\mathrm{L}})$. For the
parameters in Fig. 11, bistability occurs for $T_{\mathrm{L}}\lesssim 7.5$ K.

From these considerations, it will be relevant to derive the evolution of the
emission linewidth from the first-order two-time correlation function to
understand the system better and eventually predict the observation of a
lasing-like emission spectrum. Such calculations, however, are outside the
scope of the present article and are left for a later study.

##  IV Conclusions

We have developed a model to describe the interaction between light and matter
in a 2D material system coupled with a nanoscale electromagnetic resonator.
Light-matter interaction breaks the translational symmetry of excitons in the
2D sheet, leading to the emergence of a localized polariton state. Our study
reveals that the formation of a macroscopic population in the localized
polariton state is influenced by correction terms addressing the breakdown of
translational symmetry. As the spatial width of the cavity mode increases,
Coulomb scattering rates rise, while the correction terms saturate,
progressively exerting less influence on the system’s dynamics. Numerical
simulations demonstrate the emergence of a threshold, bistability, and
hysteresis. These phenomena can be tuned by the initial lattice temperature
while keeping the in-plane cavity size and light-matter interaction strength
constant. It is noteworthy that we observe hysteresis under non-resonant
pumping, which is distinct from the hysteresis observed in larger polariton
microcavities under resonant pumping. Thus, the cause for bistability in this
case is attributed to different factors.

##  V ACKNOWLEDGMENTS

The authors thank George Kountouris for providing Fig. 2. This paper was
supported by the Danish National Research Foundation through NanoPhoton-
Centerthe for Nanophotonics, Grant No. DNRF147.
E. V. D. acknowledges support from Independent Research Fund Denmark through
an International Postdoc Fellowship (Grant No. 0164-00014B)

##  Appendix A Expectation values

The scattering rates are computed from Wick’s theorem by expanding the six-
operator expectation value in Eq.(11). For this computation, we make the
following remarks: a) the 2D sheet is sufficiently big such that each exciton
reaction coordinate coefficient can be mapped into a momentum $\mathbf{k}$. It
follows that $U_{i,\alpha,\mathbf{k}}=\delta_{i,\mathbf{k}}$. b) The Coulomb
interaction matrix is assumed to be real
$W^{*}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}=W_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}$.
c) The approximation $W_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}\simeq
W_{000}$ still holds. d) The matrices for the change of basis in the ERC are
unitary, giving a delta function that allows us to consider only some of the
terms in the long expansion of Wick’s theorem. e) Effects of intervalley
Coulomb exchange are neglected [19]. We can then expand

$\begin{split}\footnotesize&\left<F^{\dagger}_{0}F^{\prime}_{0}\right>=\sum_{ijhi^{\prime}j^{\prime}h^{\prime}}\mathcal{U}_{i,j,h}\mathcal{U}^{*}_{i^{\prime},j^{\prime},h^{\prime}}\left<B^{\dagger}_{i}B_{j}B_{h}B^{\prime\dagger}_{h^{\prime}}B^{\prime\dagger}_{j^{\prime}}B^{\prime}_{i^{\prime}}\rho_{\mathrm{R}}\right>=\\\
&=\hbar^{2}\sum_{\begin{subarray}{c}\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}\\\
\mathbf{k}_{4},\mathbf{q},\mathbf{q}^{\prime}\\\
ijhi^{\prime}j^{\prime}h^{\prime}\end{subarray}}e^{\frac{i(t-s)}{\hbar}(\epsilon_{h^{\prime}}+\epsilon_{j^{\prime}}-\epsilon_{i^{\prime}})}W^{*}_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{q}}W_{\mathbf{k}_{3},\mathbf{k}_{4},\mathbf{q}^{\prime}}\\\
&U_{0,\mathbf{k}_{1}+\mathbf{q}}U_{i,\mathbf{k}_{2}-\mathbf{q}}U^{*}_{j,\mathbf{k}_{1}}U^{*}_{h,\mathbf{k}_{2}}U_{h^{\prime},\mathbf{k}_{3}+\mathbf{q}}U_{j^{\prime},\mathbf{k}_{4}-\mathbf{q}}U^{*}_{i^{\prime},\mathbf{k}_{3}}U^{*}_{0,\mathbf{k}_{4}}\\\
&\big{[}\delta_{ij}N_{i}\big{[}(1+N_{h})N_{j}^{\prime}\delta_{hh^{\prime}}\delta_{j^{\prime}i^{\prime}}+(1+N_{h})N_{h}^{\prime}\delta_{hj^{\prime}}\delta_{h^{\prime}i^{\prime}}\big{]}+\\\
&+\delta_{ih}N_{i}\big{[}(1+N_{j})N_{j^{\prime}}\delta_{jh^{\prime}}\delta_{j^{\prime}i^{\prime}}+(1+N_{j})N_{h^{\prime}}\delta_{jj^{\prime}}\delta_{h^{\prime}i^{\prime}}\big{]}+\\\
&+\delta_{ii^{\prime}}N_{i}\big{[}(1+N_{j})(1+N_{h})\delta_{jh^{\prime}}\delta_{j^{\prime}h}+\\\
&+(1+N_{j})(1+N_{h})\delta_{jj^{\prime}}\delta_{hh^{\prime}}\big{]}\big{]}\end{split}$
(32)

The same expansion is performed for the term,
$\left<F_{0}F^{{}^{\prime}\dagger}_{0}\right>$. It follows that the expression
for the rate of injection into the lower polariton state is

$\begin{split}\footnotesize&\left<F_{0}F^{\dagger}_{0}\right>\Big{(}e^{i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}}+e^{-i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}}\Big{)}=\\\
&(e^{i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}}+e^{-i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}})\hbar^{2}\sum_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}}\big{[}\big{[}N_{\mathbf{k}_{1}}N_{\mathbf{k}_{2}}N_{\mathbf{k}_{3}}|U_{0,\mathbf{k}_{1}}|^{2}\\\
&(W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{0}}W_{\mathbf{k}_{3},\mathbf{k}_{1},\mathbf{0}}+W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{0}}W_{\mathbf{k}_{3},\mathbf{k}_{1},\mathbf{k}_{1}-\mathbf{k}_{3}})+\\\
&+\big{[}N_{\mathbf{k}_{1}}N_{\mathbf{k}_{2}}N_{\mathbf{k}_{3}}|U_{0,\mathbf{k}_{2}}|^{2}(W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{2}-\mathbf{k}_{1}}W_{\mathbf{k}_{3},\mathbf{k}_{2},\mathbf{0}}+\\\
&+W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{2}-\mathbf{k}_{1}}W_{\mathbf{k}_{3},\mathbf{k}_{2},\mathbf{k}_{2}-\mathbf{k}_{3}})\big{]}+\big{[}N_{\mathbf{k}_{1}}N_{\mathbf{k}_{2}}(1+N_{\mathbf{k}_{3}})\\\
&|U_{0,\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{k}_{3}}|^{2}(W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{2}-\mathbf{k}_{3}}W_{\mathbf{k}_{3},\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{k}_{3},\mathbf{k}_{1}-\mathbf{k}_{3}}\\\
&+W_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{2}-\mathbf{k}_{3}}W_{\mathbf{k}_{3},\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{k}_{3},\mathbf{k}_{2}-\mathbf{k}_{3}})\big{]}\big{]}\end{split}$
(33)

A similar expression is recovered for the out-scattering rate from
$\left<F^{\dagger}_{0}F_{0}\right>\Big{(}e^{i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}}+e^{-i\frac{(t-s)\epsilon_{\mathrm{lp}}}{\hbar}}\Big{)}$.

##  Appendix B Scattering rates

We follow the steps performed in Ref. [19] to simplify the expression for
$U_{0\mathbf{k}}$ of the 2D sheet. The coupling term can be expanded into

$\begin{split}G_{0}=&\sqrt{\sum_{\mathbf{k}}{|g_{\mathbf{k}}|^{2}}}\rightarrow\\\
&\hbar
g_{\mathbf{k}}=-\frac{e_{0}}{m_{0}}\sqrt{\frac{\hbar}{2\epsilon_{0}\omega_{\mathrm{c}}}}\matrixelement{\Phi_{\mathbf{k}}}{\sum_{i}\tilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{r}_{i})\hat{\mathbf{p}}_{i}}{0}.\end{split}$
(34)

$g_{\mathbf{k}}$ is the interaction strength between an exiton and the
resonant electromagnetic field evaluated in the fermionic space. $e_{0}$ and
$m_{0}$ are the free electron charge and mass, $\epsilon_{0}$ is the electric
permittivity of free space, $\omega_{\mathrm{c}}$ is the frequency of the
cavity mode, $\Phi_{\mathbf{k}}$ is the Wannier-Mott exciton states defined as
a momentum superposition of electron-hole pairs. The mode profile for the
electric field mode of the optical cavity is
$\tilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{r}_{i},t)=\tilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{r}_{i})\delta(t)$,
which defines local and non retarded coupling dynamics, $\hat{\mathbf{p}}_{i}$
is the particle momentum. Solving the bracket in Eq. (34) through the Slater-
Condon rules [45] and making some further manipulations [19], we get to the
expression for the coupling strength

$\hbar
g_{\mathbf{k}}=-\frac{e_{0}}{m_{0}}\sqrt{\frac{\hbar}{\pi\epsilon_{0}\omega_{\mathrm{c}}a^{2}_{\mathrm{B}}S}}\int
d^{2}\mathbf{r}e^{-i\mathbf{k}\mathbf{r}}\tilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{r},z_{0})\mathbf{p}_{\mathrm{cv}}.$
(35)

$S$ is the surface area of the 2D sheet, $a_{\mathrm{B}}$ is the Bohr radius.
The 2D sheet is located at $z=z_{0}$, and the integral is performed over the
infinite extent of the 2D material. Finally,
$\mathbf{p}_{\mathrm{cv}}=\int_{V_{\mathrm{uc}}}d^{3}\mathbf{r}u^{*}_{\mathrm{c}}\hat{\mathbf{p}}u_{\mathrm{v}}(\mathbf{r})$
is the matrix element of the exciton mode, with $u_{i}(\mathbf{r})$ a Bloch
function, the subscript $i=\mathrm{c},\mathrm{v}$ the conduction and valence
band and $V_{\mathrm{uc}}$ the volume of the unit cell. From the exciton
reaction coordinates

$G_{0}=\sqrt{\sum_{\mathbf{k}}{|g_{\mathbf{k}}|^{2}}}=\sqrt{\int
J(\mathbf{\omega})d\mathbf{\omega}},$ (36)

where $J(\mathbf{\omega})$ is the spectral density of the exciton states of
the residual environment coupled to the exciton reaction coordinate $B_{0}$
[19]. When combining Eq. (35) with Eq. (36), the coupling strength becomes

$\hbar G_{0}=\sqrt{\frac{\hbar
e^{2}_{0}}{\pi\epsilon_{0}m^{2}_{0}\omega_{\mathrm{c}}a^{2}_{\mathrm{B}}}\int
d^{2}\mathbf{r}|\tilde{F}_{\mathrm{c}}(\mathbf{r},z_{0})\mathbf{p}_{\mathrm{cv}}|^{2}}.$
(37)

This expression can be widely simplified for a localized and separable mode
profile

$\hbar G_{0}=\sqrt{\frac{\hbar
e^{2}_{0}|\mathbf{n}\mathbf{p}_{\mathrm{cv}}|^{2}}{\pi\epsilon_{0}m^{2}_{0}\omega_{\mathrm{c}}a^{2}_{\mathrm{B}}L_{z}}}.$
(38)

$\mathbf{n}$ is the unit polarization vector of the mode, and the out-of-plane
confinement length

$L_{z}=\frac{\epsilon_{\mathrm{eff}}\int
d_{z}|\tilde{F}(z)|^{2}}{|\tilde{F}_{z}(z_{0})|^{2}}$ (39)

obeying the normalization requirement $\epsilon_{\mathrm{eff}}\int
d_{z}|\tilde{F}(z)|^{2}=1$, with $\epsilon_{\mathrm{eff}}$ an effective
dielectric constant to account for the dielectric response of the surrounding
material.

Also, if the later field distribution defined by the electromagnetic resonator
is assumed to be Gaussian with a confinement length scale $L$

$F_{//}(x,y)=\frac{e^{-(x^{2}+y^{2})/(2L^{2})}}{L\sqrt{\pi}},$ (40)

the exciton-mode coupling term reads

$g_{\mathbf{k}}=-\sqrt{\frac{4e^{2}_{0}|\mathbf{n}\mathbf{p}_{\mathrm{cv}}|^{2}L^{2}}{\hbar\epsilon_{0}m^{2}_{0}\omega_{\mathrm{c}}a^{2}_{\mathrm{B}}S}}\tilde{F}_{z}(z_{0})e^{-\frac{1}{2}kL}.$
(41)

It follows that the expression for
$|U_{0,\mathbf{k}}|^{2}\sim\frac{L^{2}}{S}e^{-L^{2}k^{2}}$.

From these results, converting the sum over $\mathbf{k}$ in Eq. (16) into an
integral $\sum_{\mathbf{k}}\rightarrow\frac{S}{(2\pi)^{2}}\int
d^{2}\mathbf{k}$ and solving, the following expression for the rate equations
are recovered

$\begin{split}W_{\mathrm{in}}=&\frac{4\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}\rho_{x}}{2\pi(k_{\mathrm{B}}T_{x})^{2}}\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}\mathrm{d}\theta_{2}\mathrm{d}\theta_{3}\\\
&e^{-2L^{2}(k^{\prime}_{1}k_{2}cos\theta_{2}-k^{\prime}_{1}k_{3}cos\theta_{3}}e^{-L^{2}k^{{}^{\prime}2}_{1}}e^{-\beta(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}})}\\\
&\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}),\end{split}$
(42)

Figure 13: Numerically computed out-scattering rate,
$W_{\mathrm{out}_{\mathrm{1}}}$, for $G_{0}=3.8\ \mathrm{meV}$ and different
optical cavity in-plane confinement lengths. The black dotted line is the
result expected for a quantum-well model with translational invariance, in
which case it is possible to derive an analytic expression for
$W_{\mathrm{out}_{\mathrm{1}}}$ [10]. Figure 14: Numerically computed out-
scattering rate, $W_{\mathrm{out}_{\mathrm{2}}}$, for $G_{0}=3.8\
\mathrm{meV}$ and different optical cavity in-plane confinements.

$\begin{split}W_{\mathrm{out}_{\mathrm{1}}}&=\frac{4\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}\rho^{2}_{x}}{(2\pi)^{2}k_{\mathrm{B}}T_{x}}\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}\mathrm{d}\theta_{2}\mathrm{d}\theta_{3}\\\
&e^{-L^{2}|\mathbf{k^{\prime}}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}|^{2}}e^{-\beta\epsilon_{\mathbf{k}_{3}}}\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}).\end{split}$
(43)

The variation of $W_{\mathrm{out}_{\mathrm{1}}}(T)$ with the temperature is
plotted in Fig. 13.

$\begin{split}W_{\mathrm{out}_{\mathrm{2}}}&=\frac{4\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}\rho_{x}}{2\pi(k_{\mathrm{B}}T_{x})^{2}}\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}\mathrm{d}\theta_{2}\mathrm{d}\theta_{3}\\\
&e^{-L^{2}|\mathbf{k^{\prime}}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}|^{2}}\big{(}e^{-\beta(\epsilon_{\mathrm{lp}}+2\epsilon_{\mathbf{k}_{3}}-\epsilon_{\mathbf{k}_{2}})}+e^{-\beta(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathbf{k}_{2}})}\big{)}\\\
&\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}),\end{split}$
(44)

whose evolution as a function of the temperature is plotted in Fig. 14.

$\begin{split}W_{\mathrm{in_{c}}}&=\frac{8\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}(2\pi)^{2}}{(k_{\mathrm{B}}T_{x})^{3}}\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}e^{-L^{2}k^{\prime
2}_{1}}\\\
&e^{-\beta(2\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}})}\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}),\\\
\end{split}$ (45)
$\begin{split}W_{\mathrm{out_{c}}}&=\frac{8\hbar|W_{000}|^{2}|X_{\mathrm{lp}}|^{2}L^{2}S^{2}2\pi\rho_{x}}{(k_{\mathrm{B}}T_{x})^{2}}\int\mathrm{d}\epsilon_{\mathbf{k}_{2}}\mathrm{d}\epsilon_{\mathbf{k}_{3}}e^{-L^{2}k^{{}^{\prime}2}_{1}}\\\
&e^{-\beta(\epsilon_{\mathbf{k}_{2}}+\epsilon_{\mathbf{k}_{3}})}\mathcal{H}(\epsilon_{\mathbf{k}_{3}}+\epsilon_{\mathrm{lp}}-\epsilon_{\mathbf{k}_{2}}),\end{split}$
(46)

with

$\begin{split}e&{}^{-L^{2}|\mathbf{k^{\prime}}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}|^{2}}=e^{-L^{2}(k^{2}_{2}+k^{2}_{3}-2k_{2}k_{3}cos(\theta_{3}-\theta_{2}))}\\\
&\ \ e^{-L^{2}(k^{\prime
2}_{1}+2k^{\prime}_{1}k_{2}cos\theta_{2}-2k^{\prime}_{1}k_{3}cos\theta_{3}})\\\
k&{}^{\prime
2}=k^{2}_{3}-k^{2}_{2}+\frac{\epsilon_{\mathrm{lp}}2m}{\hbar^{2}}\end{split}$
(47)

and $\mathcal{H}$ the Heaviside function. The third-order term from the out-
scattering rate $\sim N^{3}_{x}N_{0}$ can be neglected here. The evolution of
$W_{\mathrm{out}_{\mathrm{c}}}(T)$ as a function of the temperature is plotted
in Fig. 15.

Figure 15: Numerically computed out-scattering rate,
$W_{\mathrm{out}_{\mathrm{c}}}$, for $G_{0}=3.8\ \mathrm{meV}$ and different
optical cavity in-plane confinements.

Comparing the temperature-dependence of the scattering rates shown in Fig. 5,
7, 13, 14, 15, we observe that the maximum value for the scattering rate is
reached at a lower temperature for the in-scattering rates compared to the
out-scattering rates. The letters keep increasing with $T_{x}$. The Maxwell-
Boltzmann distribution broadens at higher temperatures, and particles have
more thermal energy. This, in turn, means that particles are more likely to
have higher kinetic energies and occupy higher energy states. This increased
energy makes it easier for the particles to scatter out of the lower polariton
state.

## References

  * Li  _et al._ [2020] Y. Li, H. Sun, L. Gan, J. Zhang, J. Feng, D. Zhang, and C.-Z. Ning, Optical properties and light-emission device applications of 2-D layered semiconductors, Proc. IEEE Inst. Electr. Electron. Eng. 108, 676 (2020).
  * Mueller and Malic [2018] T. Mueller and E. Malic, Exciton physics and device application of two-dimensional transition metal dichalcogenide semiconductors, npj 2D Materials and Applications 2, 10.1038/s41699-018-0074-2 (2018).
  * Deng  _et al._ [2002] H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Yamamoto, Condensation of semiconductor microcavity exciton polaritons, Science 298, 199–202 (2002).
  * Kasprzak  _et al._ [2006] 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 L. S. Dang, Bose–einstein condensation of exciton polaritons, Nature 443, 409 (2006).
  * Deng  _et al._ [2007] H. Deng, G. S. Solomon, R. Hey, K. H. Ploog, and Y. Yamamoto, Spatial coherence of a polariton condensate, Physical Review Letters 99, 10.1103/physrevlett.99.126403 (2007).
  * Anton-Solanas  _et al._ [2021] C. Anton-Solanas, M. Waldherr, M. Klaas, H. Suchomel, T. H. Harder, H. Cai, E. Sedov, S. Klembt, A. V. Kavokin, S. Tongay, K. Watanabe, T. Taniguchi, S. Höfling, and C. Schneider, Bosonic condensation of exciton polaritons in an atomically thin crystal, Nature Materials 20, 1233 (2021).
  * Waldherr  _et al._ [2018] M. Waldherr, N. Lundt, M. Klaas, S. Betzold, M. Wurdack, V. Baumann, E. Estrecho, A. Nalitov, E. Cherotchenko, H. Cai, E. A. Ostrovskaya, A. V. Kavokin, S. Tongay, S. Klembt, S. Höfling, and C. Schneider, Observation of bosonic condensation in a hybrid monolayer MoSe2-GaAs microcavity, Nature Communications 9, 10.1038/s41467-018-05532-7 (2018).
  * Scafirimuto  _et al._ [2017] F. Scafirimuto, D. Urbonas, U. Scherf, R. F. Mahrt, and T. Stöferle, Room-temperature exciton-polariton condensation in a tunable zero-dimensional microcavity, ACS Photonics 5, 85 (2017).
  * Tassone and Yamamoto [1999] F. Tassone and Y. Yamamoto, Exciton-exciton scattering dynamics in a semiconductor microcavity and stimulated scattering into polaritons, Physical Review B 59, 10830 (1999).
  * Porras  _et al._ [2002] D. Porras, C. Ciuti, J. J. Baumberg, and C. Tejedor, Polariton dynamics and bose-einstein condensation in semiconductor microcavities, Physical Review B 66, 10.1103/physrevb.66.085304 (2002).
  * Tartakovskii  _et al._ [2000] A. I. Tartakovskii, M. Emam-Ismail, R. M. Stevenson, M. S. Skolnick, V. N. Astratov, D. M. Whittaker, J. J. Baumberg, and J. S. Roberts, Relaxation bottleneck and its suppression in semiconductor microcavities, Physical Review B 62, R2283 (2000).
  * Kountouris  _et al._ [2022] G. Kountouris, J. Mørk, E. Denning, and P. Kristensen, Modal properties of dielectric bowtie cavities with deep sub-wavelength confinement, Optics Express 30, 40367 (2022).
  * Albrechtsen  _et al._ [2022] M. Albrechtsen, B. V. Lahijani, R. E. Christiansen, V. T. H. Nguyen, L. N. Casses, S. E. Hansen, N. Stenger, O. Sigmund, H. Jansen, J. Mørk, and S. Stobbe, Nanometer-scale photon confinement in topology-optimized dielectric cavities, Nature Communications 13, 10.1038/s41467-022-33874-w (2022).
  * Garg  _et al._ [1985] A. Garg, J. N. Onuchic, and V. Ambegaokar, Effect of friction on electron transfer in biomolecules, J. Chem. Phys. 83, 4491 (1985).
  * Thoss  _et al._ [2001] M. Thoss, H. Wang, and W. H. Miller, Self-consistent hybrid approach for complex systems: Application to the spin-boson model with debye spectral density, J. Chem. Phys. 115, 2991 (2001).
  * Iles-Smith  _et al._ [2014] J. Iles-Smith, N. Lambert, and A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems, Physical Review A 90, 10.1103/physreva.90.032114 (2014).
  * Imamoglu  _et al._ [1996] A. Imamoglu, R. J. Ram, S. Pau, and Y. Yamamoto, Nonequilibrium condensates and lasers without inversion: Exciton-polariton lasers, Physical Review A 53, 4250 (1996).
  * Snoke [2012] D. Snoke, Polariton condensation and lasing, in  _Springer Series in Solid-State Sciences_ (Springer Berlin Heidelberg, 2012) p. 307–327.
  * Denning  _et al._ [2022a] E. V. Denning, M. Wubs, N. Stenger, J. Mørk, and P. T. Kristensen, Quantum theory of two-dimensional materials coupled to electromagnetic resonators, Physical Review B 105, 10.1103/physrevb.105.085306 (2022a).
  * Valatin [1958] J. G. Valatin, Comments on the theory of superconductivity, Il Nuovo Cimento 7, 843 (1958).
  * Bogoljubov  _et al._ [1958] N. N. Bogoljubov, V. V. Tolmachov, and D. V. Širkov, On a new method in the theory of superconductivity, Il Nuovo Cimento 7, 794 (1958).
  * Ching  _et al._ [1998] E. S. C. Ching, P. T. Leung, A. M. van den Brink, W. M. Suen, S. S. Tong, and K. Young, Quasinormal-mode expansion for waves in open systems, Reviews of Modern Physics 70, 1545 (1998).
  * Kristensen and Hughes [2013] P. T. Kristensen and S. Hughes, Modes and mode volumes of leaky optical cavities and plasmonic nanoresonators, ACS Photonics 1, 2 (2013).
  * Kristensen  _et al._ [2020] P. T. Kristensen, K. Herrmann, F. Intravaia, and K. Busch, Modeling electromagnetic resonators using quasinormal modes, Advances in Optics and Photonics 12, 612 (2020).
  * Franke  _et al._ [2019] S. Franke, S. Hughes, M. K. Dezfouli, P. T. Kristensen, K. Busch, A. Knorr, and M. Richter, Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics, Physical Review Letters 122, 10.1103/physrevlett.122.213901 (2019).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of MoS2 and other group-VI dichalcogenides, Physical Review Letters 108, 10.1103/physrevlett.108.196802 (2012).
  * Cao  _et al._ [2012] T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, Q. Niu, P. Tan, E. Wang, B. Liu, and J. Feng, Valley-selective circular dichroism of monolayer molybdenum disulphide, Nature Communications 3, 10.1038/ncomms1882 (2012).
  * Li  _et al._ [2014] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Measurement of the optical dielectric function of monolayer transition-metal dichalcogenides: Mos2, mose2, ws2, and wse2, Physical Review B 90, 10.1103/physrevb.90.205422 (2014).
  * Schneider  _et al._ [2018] C. Schneider, M. M. Glazov, T. Korn, S. Höfling, and B. Urbaszek, Two-dimensional semiconductors in the regime of strong light-matter coupling, Nature Communications 9, 10.1038/s41467-018-04866-6 (2018).
  * Lundt  _et al._ [2016] N. Lundt, A. Maryński, E. Cherotchenko, A. Pant, X. Fan, S. Tongay, G. Sek, A. V. Kavokin, S. Höfling, and C. Schneider, Monolayered MoSe2: a candidate for room temperature polaritonics, 2D Materials 4, 015006 (2016).
  * Al-Ani  _et al._ [2021] I. A. M. Al-Ani, K. As'Ham, Lujun, A. E. Miroshnichenko, and H. T. Hattori, Enhanced strong coupling of TMDC monolayers by bound state in the continuum, Laser & Photonics Reviews 15, 2100240 (2021).
  * Rochat  _et al._ [2000] G. Rochat, C. Ciuti, V. Savona, C. Piermarocchi, A. Quattropani, and P. Schwendimann, Excitonic bloch equations for a two-dimensional system of interacting excitons, Physical Review B 61, 13856 (2000).
  * Verger  _et al._ [2006] A. Verger, C. Ciuti, and I. Carusotto, Polariton quantum blockade in a photonic dot, Physical Review B 73, 10.1103/physrevb.73.193306 (2006).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (OXFORD UNIVERSITY PRESS, 2002).
  * Porras and Tejedor [2003] D. Porras and C. Tejedor, Linewidth of a polariton laser: theoretical analysis of self-interaction effects, Physical Review B 67, 10.1103/physrevb.67.161310 (2003).
  * Wick [1950] G. C. Wick, The evaluation of the collision matrix, Physical Review 80, 268 (1950).
  * Deng  _et al._ [2010] H. Deng, H. Haug, and Y. Yamamoto, Exciton-polariton bose-einstein condensation, Reviews of Modern Physics 82, 1489 (2010).
  * Note [1] The $i_{th}$ exciton state occupation number is $N^{x}_{i}=N^{x}\exp{-\beta\epsilon_{i}}$, with $\beta=(T_{x}k_{\mathrm{B}})^{-1}$ and $\epsilon_{i}$ the corresponding energy.
  * Sivakumar [2018] S. Sivakumar, Microcavities, by alexey v. kavokin, jeremy j. baumberg, guillaume maupuech, and fabrice p. laussy, Contemp. Phys. 59, 203 (2018).
  * Denning  _et al._ [2022b] E. V. Denning, M. Wubs, N. Stenger, J. Mørk, and P. T. Kristensen, Cavity-induced exciton localization and polariton blockade in two-dimensional semiconductors coupled to an electromagnetic resonator, Physical Review Research 4, 10.1103/physrevresearch.4.l012020 (2022b).
  * Selig  _et al._ [2016] M. Selig, G. Berghäuser, A. Raja, P. Nagler, C. Schüller, T. F. Heinz, T. Korn, A. Chernikov, E. Malic, and A. Knorr, Excitonic linewidth and coherence lifetime in monolayer transition metal dichalcogenides, Nature Communications 7, 10.1038/ncomms13279 (2016).
  * Huang  _et al._ [2000] R. Huang, F. Tassone, and Y. Yamamoto, Experimental evidence of stimulated scattering of excitons into microcavity polaritons, Physical Review B 61, R7854 (2000).
  * Baas  _et al._ [2004] A. Baas, J.-P. Karr, M. Romanelli, A. Bramati, and E. Giacobino, Optical bistability in semiconductor microcavities in the nondegenerate parametric oscillation regime: Analogy with the optical parametric oscillator, Phys. Rev. B Condens. Matter Mater. Phys. 70 (2004).
  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen, The road towards polaritonic devices, Nature Materials 15, 1061 (2016).
  * Grosso and Parravicini [2013] G. Grosso and G. Parravicini,  _Solid State Physics_ (Academic Press, Oxford, UK, 2013).