36 |  {Find(Is(Object(guitar)), FaceObject, GetAbove) $\rightarrow$ Sharpen} |  Objects |  Sharpen faces of people playing guitar. |  6
37 |  {Find(Is(Word(319)), Object(car),
GetParents) $\rightarrow$ Blur} |  Objects |  Blur car with number 319. |  7
38 |  {Union(Is(Object(car)), Is(Object(bicycle))) $\rightarrow$ Brighten} |  Objects |  Brighten all cars and bicycles. |  7
39 |  {Find(Is(Object(person)),
Object(bicycle), GetBelow) $\rightarrow$ Brighten} |  Objects |  Brighten all bicycles that are being ridden. |  7
40 |  {Find(Is(Object(bicycle)),
BelowAge(18), GetAbove) $\rightarrow$ Blur} |  Objects |  Blur the faces of children riding bicycles. |  7
41 |  {Complement(Union(Is(Object(car)), Is(Object(bicycle)))) $\rightarrow$ Blackout} |  Objects |  Blackout all objects except cars and bicycles. |  8
42 |  {Intersect(Is(TextObject),
Complement(Filter(Is(Object(car)),
TextObject))) $\rightarrow$ Blackout} |  Objects |  Blackout all text not on a car. |  9
43 |  {Union(Is(Object(bicycle)), Is(Object(car)), Is(Object(person))) $\rightarrow$ Brighten} |  Objects |  Brighten all bicycles, cars, and people. |  10
44 |  {Intersect(Is(FaceObject),
Complement(Find(Is(Object(bicycle)), FaceObject, GetAbove))) $\rightarrow$ Blur} |  Objects |  Blur faces of people not riding bicycles. |  10
45 |  {Union(Is(Object(guitar)), Find(Is(Object(guitar)), FaceObject, GetAbove)) $\rightarrow\textsf{Brighten}$} |  Objects |  Brighten all guitars and people playing guitar. |  10
46 |  {Intersect(Is(Face),
Complement(Find(Is(Object(guitar)), FaceObject, GetAbove))) $\rightarrow$ Blur} |  Objects |  Blur faces of people not playing guitar. |  10
47 |  {Intersect(Is(Object(bicycle)),
Complement(Find(Is(Object(person)), Object(bicycle), GetBelow))) $\rightarrow$ Sharpen} |  Objects |  Sharpen bicycles that are not being ridden. |  12
48 |  {Intersect(Is(Object(bicycle)), Complement(Find(Is(BelowAge(18)), Object(bicycle), GetBelow))) $\rightarrow$ Sharpen} |  Objects |  Sharpen all bicycles that are not ridden by a child. |  12
49 |  {Intersect(Is(Object(cat)),
Complement(Find(Is(Object(cat)),
Object(cat), GetBelow))) $\rightarrow$ Crop } |  Objects |  Crop image to feature just topmost cat. |  12
50 |  {Intersect(Find(Is(Object(cat)),
Object(cat), GetRight),
Find(Is(Object(cat)), Object(cat),
GetLeft)) $\rightarrow$ Brighten} |  Objects |  Brighten cats that are between two other cats. |  15

##  Appendix C List of Neural Attributes

  * •

FaceObject. Objects identified by Amazon Rekognition’s face recognition model.

  * •

Smiling. Faces identified as smiling by Rekognition’s face recognition model.

  * •

EyesOpen. Faces identified as having open eyes by Rekognition’s face
recognition model.

  * •

MouthOpen. Face identified as having open mouths by Rekognition’s face
recognition model.

  * •

BelowAge$(N)$. Rekognition’s face recognition model returns an upper and lower
bound on a face’s age. This predicate matches a face if the upper bound on its
age is less than $N$.

  * •

Face$(N)$. Rekognition assigns each unique face an identified. This predicate
matches a face if its identified is equal to $N$.

  * •

TextObject. Objects identified by Amazon Rekognition’s text recognition model.

  * •

Word$(W)$. Rekognition’s text recognition model returns the text value of each
identified text object. This predicate matches a text object if its text value
equals $W$.

  * •

PhoneNumber. Text objects whose value matches the format of a North American
phone number.

  * •

Price. Text objects whose value matches the format of a price.

  * •

Object$(O)$. Objects identified by Rekognition’s object recognition model.
This model identifies 238 unique types of objects. This predicate matches an
object whose type is equal to $O$.


# Exact analytical solution for density matrix of a non-equilibrium polariton
Bose-Einstein condensate

Vladislav Yu. Shishkov<EMAIL_ADDRESS>Evgeny S. Andrianov  Dukhov
Research Institute of Automatics (VNIIA), 22 Sushchevskaya, Moscow 127055,
Russia  Moscow Institute of Physics and Technology, 9 Institutskiy pereulok,
Dolgoprudny 141700, Moscow region, Russia  Skolkovo Institute of Science and
Technology, Moscow, Russian Federation  Anton V. Zasedatelev  Skolkovo
Institute of Science and Technology, Moscow, Russian Federation  Pavlos G.
Lagoudakis  Skolkovo Institute of Science and Technology, Moscow, Russian
Federation  Department of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom  Yurii E. Lozovik  Institute for
Spectroscopy RAS, 5 Fizicheskaya, Troitsk 142190, Russia  Moscow Institute of
Electronics and Mathematics, National Research University Higher School of
Economics, 101000, Moscow, Russia  Skolkovo Institute of Science and
Technology, Moscow, Russian Federation

###### Abstract

In this letter, we give an analytical quantum description of a non-equilibrium
polariton Bose-Einstein condensate (BEC) based on the solution of the master
equation for the full polariton density matrix in the limit of fast
thermalization. We find the density matrix of a non-equilibrium BEC, that
takes into account quantum correlations between all polariton states. We show
that the formation of BEC is accompanied by the build-up of cross-correlations
between the ground state and the excited states reaching their highest values
at the condensation threshold. Despite the non-equilibrium nature of polariton
systems, we show the average population of polariton states exhibits the Bose-
Einstein distribution with an almost zero effective chemical potential above
the condensation threshold similar to an equilibrium BEC. We demonstrate that
above threshold the effective temperature of polariton condensate drops below
the reservoir temperature.

††preprint: APS/123-QED

Introduction. The theoretical prediction of BEC was made at the dawn of modern
statistical physics [1, 2]. However, the first experimental realization of the
BEC required more than 70 years of physicists’ efforts [3, 4]. BEC was first
observed for atoms in a magneto-optical trap at cryogenic temperatures [3, 4].
A rapid development of semiconductor technologies made it possible to
experimentally realize polariton BEC in solid state at elevated temperature
[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 12, 15, 16, 17, 18, 19]. Polaritons are
hybrid quasiparticles formed by the strong light-matter interaction of photons
in a cavity with excitons of a semiconductor material or organic dyes [7, 20].
Due to their small effective mass and low density of states [8, 9, 14, 15],
the critical temperature for a polariton BEC can be as high as tens or even
hundreds of kelvin [5, 6, 19]. Such light-matter condensed states attract
fundamental as well as practical interest due to their potential in all-
optical manipulation and low-energy optoelectronic applications [21, 16].

Being open dissipative systems polariton condensates intrinsically are far
from equilibrium. To maintain the number of polaritons at a certain level an
external pumping is required. Therefore, the methods of equilibrium
thermodynamics must be applied with caution to the description of non-
equilibrium polariton condensates [7, 22, 23, 24]. The description of the
dynamics of polariton condensates require more sophisticated approaches (for a
review, see the Supplement in [16]). Generally, complete information about the
polariton condensate can be accessed through its full density matrix. When
dissipation processes are described in the Markov approximation, their
dynamics obeys the corresponding Lindblad equation [23]. However, due to the
large number of excited states, the master equation for the polariton density
matrix cannot be solved exactly. Previous studies largely rely on the
semiclassical Maxwell–Boltzmann equations to the average population of
polaritons in condensates [22, 25, 26, 27, 28] or other mean-field theories
[29, 30, 31]. This approach takes into account both dissipation and pumping as
well as the processes leading to the thermalization of the polaritons [23,
32]. The Maxwell–Boltzmann equations reproduce the observed Gibbs distribution
of the polaritons in the excited states [32, 33]. However, this approach does
not provide any information neither on the coherence of polaritons nor on the
correlations between the polariton states [23, 34]. Recently, a theory beyond
the mean-field approximation was developed [35]. This theory consistently
describes lasing in weakly and strongly coupled organic systems with strong
exciton-vibron coupling. It takes two-time correlations into account allowing
to study evolution of coherence and thermalization in such systems, albeit
limited to the first-order correlations.

A way to get information about the coherence of the polariton BEC is to reduce
the complete master equation to an effective one for the density matrix of the
ground polariton state only, where the contribution of the excited polariton
states is taken into account as a whole effectively [27, 34, 30, 36, 37, 38,
39, 40]. Such reduced equation for the density matrix allows one to obtain
complete information about the quantum state of the ground polariton state
upon the BEC formation [23, 34]. However, this approach does not provide any
information about the coherence of the polaritons in the excited states, their
statistical properties and their cross-correlations between the ground state.

In this letter, we propose a complete quantum-mechanical approach for
obtaining the density matrix of non-equilibrium polariton condensates in the
limit of fast polariton thermalization. We find the stationary density matrix
of the polariton BEC, which allows us to analyze the correlations between all
polariton states. We show that cross-correlations between polaritons in the
ground state and those in the excited states arises at the condensation
threshold. This property is a unique feature of non-equilibrium BEC. We show
that the number of polaritons in the ground and excited states obeys
Bose–Einstein distribution with effective temperature well below the reservoir
temperature and an almost zero effective chemical potential above condensation
threshold.

Master equation for the lower polaritons. Strong light-matter interaction of
the cavity photons with an optical transition of active material gives rise
for new eigenstates, namely lower and upper polariton branches. The polaritons
of the upper branch and excitons [18, 8]) are able to scatter towards lower
polariton branch forming BEC at high enough excitation density. There are
several mechanisms allowing such energy relaxation in polariton systems.
Scattering with electrons [41] acoustic [23] and high-energy optical phonons
[42] or vibrons in organic materials [18, 43] are among them. Pair-particle
scattering processes can also contribute to energy relaxation towards the low-
lying polariton states at high polariton densities [44].

At the same time polaritons undergo multiple scattering at the lower polariton
branch followed by thermalization and condensation above critical polariton
density. Since the non-equilibrium BEC occurs in the lower polariton branch,
below we consider the dynamics of the lower polariton branch only. We describe
the energy transfer from the upper polaritons or excitons to the lower
polaritons as an effective incoherent pumping which is reasonable assumption
for majority of experimental realizations [6, 19, 13]

We suppose that the polaritons with frequencies ${\omega_{j}}$ of the lower
polariton branch near the minimum are described by bosonic creation
$\hat{a}_{j}^{\dagger}$ and annihilation $\hat{a}_{j}$ operators [23, 34]
which obey the commutation relation
$\left[{{{\hat{a}}_{j}},\hat{a}_{j^{\prime}}^{\dagger}}\right]={\delta_{jj^{\prime}}}$.
In this case, the Hamiltonian of the polaritons takes the form

${\hat{H}_{{\rm{LP}}}}=\sum\limits_{j=0}^{M}{{\omega_{j}}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}}$
(1)

The state with $j=0$ is the ground state.

We describe the dynamics of the polaritons in the lower polariton branch
through the density matrix $\hat{\rho}$. We consider the relaxation and
pumping processes in the Born–Markov approximation [45]. In such a case, the
density matrix is governed by the master equation in the Lindblad form. The
dissipation of the lower polaritons is described by the Lindblad superoperator

${L_{{\rm{diss}}}}\left({\hat{\rho}}\right)=\sum\limits_{j=0}^{M}{{\gamma_{j}}\left({{{\hat{a}}_{j}}\hat{\rho}\hat{a}_{j}^{\dagger}-\frac{1}{2}\hat{\rho}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}-\frac{1}{2}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}\hat{\rho}}\right)}$
(2)

where ${\gamma_{j}}$ is the dissipation rate of the $j$-th state. The
effective incoherent pumping of the lower polaritons can be described by the
Lindblad superoperator

${L_{{\rm{pump}}}}\left({\hat{\rho}}\right)=\sum\limits_{j=0}^{M}{{\kappa_{j}}\left({{{\hat{a}}_{j}}\hat{\rho}\hat{a}_{j}^{\dagger}-\frac{1}{2}\hat{\rho}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}-\frac{1}{2}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}\hat{\rho}}\right)}\\\
+\sum\limits_{j=0}^{M}{{\kappa_{j}}\left({\hat{a}_{j}^{\dagger}\hat{\rho}{{\hat{a}}_{j}}-\frac{1}{2}\hat{\rho}{{\hat{a}}_{j}}\hat{a}_{j}^{\dagger}-\frac{1}{2}{{\hat{a}}_{j}}\hat{a}_{j}^{\dagger}\hat{\rho}}\right)}$
(3)

where ${\kappa_{j}}$ is the pumping rate of the $j$-th state. The operator
${L_{{\rm{pump}}}}\left({\hat{\rho}}\right)$ leads to the following dynamics
of the average polariton population in the $j$-th state
$\left\langle{\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}}\right\rangle$ and the
average polariton amplitude in the $j$-th state
$\left\langle{{{\hat{a}}_{j}}}\right\rangle$:
${\left({{{d\left\langle{\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}}\right\rangle}\mathord{\left/{\vphantom{{d\left\langle{\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}}\right\rangle}{dt}}}\right.\kern-1.2pt}{dt}}}\right)_{{\rm{pump}}}}={\kappa_{j}}$
and
${\left({{{d\left\langle{{{\hat{a}}_{j}}}\right\rangle}\mathord{\left/{\vphantom{{d\left\langle{{{\hat{a}}_{j}}}\right\rangle}{dt}}}\right.\kern-1.2pt}{dt}}}\right)_{{\rm{pump}}}}=0$.
Therefore, the action of incoherent pumping (3) leads to the excitation of a
certain number of polaritons per unit time in the corresponding state, without
affecting their phase [46].

Thermalization of the lower polaritons may occur due to different physical
processes depending on the system. For example, in organic polariton systems
thermalization occurs due to their nonlinear interaction with low frequency
vibrations [31, 47, 48, 49, 43, 50, 51, 52]. For polariton states in inorganic
semiconductors thermalization predominantly goes through interactions with
acoustic phonons or free charges [23]. Also the thermalization of the
polaritons can occur due to polariton-polariton scattering [14]. However,
regardless of the mechanism, this thermalization can be described through the
Lindblad superoperator [23, 34].

${L_{{\rm{therm}}}}\left({\hat{\rho}}\right)=\\\
\sum\limits_{j=0}^{M}{\sum\limits_{k=0}^{M}{{\Gamma_{jk}}\left({{{\hat{a}}_{j}}\hat{a}_{k}^{\dagger}\hat{\rho}{{\hat{a}}_{k}}\hat{a}_{j}^{\dagger}-\frac{1}{2}\hat{\rho}{{\hat{a}}_{k}}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}\hat{a}_{k}^{\dagger}-\frac{1}{2}{{\hat{a}}_{k}}\hat{a}_{j}^{\dagger}{{\hat{a}}_{j}}\hat{a}_{k}^{\dagger}\hat{\rho}}\right)}}$
(4)

where ${\Gamma_{jk}}$ is the transition rate from the $j$-th polariton state
to the $k$-th state. The thermalization rates ${\Gamma_{jk}}$ obey the
Kubo–Martin–Schwinger relation
${\Gamma_{jk}}/{\Gamma_{kj}}=\exp\left({\left({{\omega_{j}}-{\omega_{k}}}\right)/T}\right)$,
where $T$ is the temperature of intermolecular oscillations of the organic
dyes or the temperature of the phonons in the semiconductors.

Thus, the complete master equation for the density matrix of the polaritons at
the lower polariton branch $\hat{\rho}$ has the form

${\frac{\partial\hat{\rho}}{\partial
t}}=\frac{i}{\hbar}\left[{\hat{\rho},{{\hat{H}}_{{\rm{LP}}}}}\right]+{L_{{\rm{diss}}}}\left({\hat{\rho}}\right)+{L_{{\rm{pump}}}}\left({\hat{\rho}}\right)+{L_{{\rm{therm}}}}\left({\hat{\rho}}\right)$
(5)

Usually it is difficult to solve the master equation (5). For instance, if we
take into account $M+1$ polariton states and $N$ excitations in each state,
then the total number of differential equations given by (5) is
${\left({N+1}\right)^{2\left({M+1}\right)}}$. Here we consider the fast
thermalization limit and show that in this case the master equation (5) can be
reduced to $\left({N+1}\right)$ differential equations.

Below we describe our approach in two stages. First, we find all possible
thermalized density matrices of the polaritons. We assume that thermalization
is the fastest process in the system, i.e.,
${\Gamma_{0j}}(1+\langle\hat{a}^{\dagger}_{0}\hat{a}_{0}\rangle)\gg{\kappa_{j}},{\gamma_{j}}$
(see Discussion). In such a case, at times
$t\ll\gamma_{j}^{-1},\kappa_{j}^{-1}$ the density matrix obeys the approximate
differential equation

$\frac{{d\hat{\rho}}}{{dt}}=\frac{i}{\hbar}\left[{\hat{\rho},{{\hat{H}}_{{\rm{LP}}}}}\right]+{L_{{\rm{therm}}}}\left({\hat{\rho}}\right)$
(6)

After the time
${\Gamma_{0j}}^{-1}(1+\langle\hat{a}^{\dagger}_{0}\hat{a}_{0}\rangle)^{-1}\ll
t\ll\gamma_{j}^{-1},\kappa_{j}^{-1}$, the system reach its quasistationary
state which is determined by the following equation:

$\frac{i}{\hbar}\left[{\hat{\rho},{{\hat{H}}_{{\rm{LP}}}}}\right]+{L_{{\rm{therm}}}}\left({\hat{\rho}}\right)=0$
(7)

In the limits
${\kappa_{j}}/{\Gamma_{0j}}(1+\langle\hat{a}^{\dagger}_{0}\hat{a}_{0}\rangle)\to
0$ and
${\gamma_{j}}/{\Gamma_{0j}}(1+\langle\hat{a}^{\dagger}_{0}\hat{a}_{0}\rangle)\to
0$, we can suppose that the system reaches quasistationary state at
instantaneous time that is defined by the condition (7). Applying the theory
[53], we can obtain the general form of the density matrix of this
quasistationary state. Indeed, the relaxation operator
${L_{{\rm{therm}}}}\left({\hat{\rho}}\right)$, defined by the expression (4),
conserves the total number of lower polaritons [23]. Therefore the operator of
total number of lower polaritons
$\sum_{j=0}^{M}\hat{a}^{\dagger}_{j}\hat{a}_{j}$ is the integral of motion for
the thermalization process. According to [53], the presence of this integral
of motion implies that the system has invariant subspaces
$\left|{{n_{0}},{n_{1}},...,{n_{M}}}\right\rangle\left\langle{{n_{0}},{n_{1}},...,{n_{M}}}\right|$
with the total number of polaritons equal to
$\sum\nolimits_{j=0}^{M}{{n_{j}}}=N$. Being in the invariant subspace at
initial moment in time, the system stays in this invariant subspace in the
subsequent moments. In each invariant subspace the Gibbs distribution is
established with a temperature $T$ [53]. The corresponding density matrix is
$Z_{N}^{-1}\sum\nolimits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}\left|{{n_{0}},...,{n_{M}}}\right\rangle\left\langle{{n_{0}},...,{n_{M}}}\right|}$,
where ${w_{j}}=\exp\left({\left({{\omega_{0}}-{\omega_{j}}}\right)/T}\right)$,
and ${Z_{N}}$ is the partition function given that the total number of
polaritons is $N$

${Z_{N}}=\sum\limits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}}$
(8)

Some properties of the partition function ${Z_{N}}$ are discussed in the
Supplemental Material I.1. The general expression for the density matrix
$\hat{\rho}\left(t\right)$ that fulfills the condition (7) is the sum of the
Gibbs distributions in each invariant subspace taken with the coefficients
${P_{N}}\left(t\right)$

$\hat{\rho}\left(t\right)=\sum\limits_{N=0}^{+\infty}{P_{N}}\left(t\right)\frac{1}{{{Z_{N}}}}\sum\limits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}}\\\
\times{\left|{{n_{0}},...,{n_{M}}}\right\rangle\left\langle{{n_{0}},...,{n_{M}}}\right|}$
(9)

The coefficients ${P_{N}}\left(t\right)$ are the probabilities that there are
$N$ polaritons in total in the low polariton branch, consequently
$\sum\nolimits_{N=0}^{+\infty}{{P_{N}}\left(t\right)}=1$.

On the second stage, we substitute the general thermalized density matrix (9)
into the master equation (5) to account for driven-dissipative behaviour of
the system. As a result, we get equations for the probabilities
${P_{N}}\left(t\right)$ (see the Supplemental Material I.2)

${\frac{\partial{P_{0}}\left(t\right)}{\partial
t}}={\frac{{d_{0}}}{{Z_{1}}}}{P_{1}}\left(t\right)-{\frac{{p_{0}}}{{Z_{0}}}}{P_{0}}\left(t\right)$
(10) ${\frac{\partial{P_{N}}\left(t\right)}{\partial
t}}={\frac{{d_{N}}}{{Z_{N+1}}}}{P_{N+1}}\left(t\right)-{\frac{\left({{d_{N-1}}+{p_{N}}}\right)}{{Z_{N}}}}{P_{N}}\left(t\right)\\\
+{\frac{{p_{N-1}}}{{Z_{N-1}}}}{P_{N-1}}\left(t\right),\,\,\,\,\text{for}\,\,N>0$
(11)

where
${d_{N}}=\sum\nolimits_{n=0}^{N}{{Z_{N-n}}\sum\nolimits_{j=0}^{M}{\left({{\gamma_{j}}+{\kappa_{j}}}\right)w_{j}^{n+1}}}$
and
${p_{N}}=\sum\nolimits_{n=0}^{N}{{Z_{N-n}}\sum\nolimits_{j=0}^{M}{{\kappa_{j}}w_{j}^{n}}}$.

In the case of steady-state non-resonant pumping the stationary solution of
Equations (10)-(11) can be written in recurrent form:

${P_{1}}={\frac{{p_{0}}}{{d_{0}}}}{\frac{{Z_{1}}}{{Z_{0}}}}{P_{0}}$ (12)
${P_{N+1}}={\frac{\left({{d_{N-1}}+{p_{N}}}\right)}{{d_{N}}}}{\frac{{Z_{N+1}}}{{Z_{N}}}}{P_{N}}-{\frac{{p_{N-1}}}{{d_{N}}}}{\frac{{Z_{N+1}}}{{Z_{N-1}}}}{P_{N-1}}$
(13)

On the basis of our theory one can find the complete density matrix of the
polaritons at the low polariton branch.

Formation of non-equilibrium BEC. In this section we apply the developed
theory to a polariton system consisting of $M+1$ states equidistant in
frequency. This condition on the frequencies of the polariton states
corresponds to a constant density of states in the continuous limit of the
two-dimensional system with quadratic dispersion. The frequency of the $j$-th
polariton state is
$\Delta{\omega_{j}}=\left({{\omega_{M}}-{\omega_{0}}}\right)\times j/M$. We
consider the decay rates to be the same for the entire lower polariton branch
and equal $\gamma$. We also assume that the incoherent pumping is
characterized by a rate $\kappa$, which does not depend on time, and acts only
on the polariton state with $j=M$.

The stationary density matrix of the polariton system under consideration has
the form (9), where the coefficients ${P_{N}}\left(t\right)$ are time-
independent given by the expressions (12)-(13). Following Eq. (12)-(13) the
stationary density matrix is well-defined by two dimensionless parameters,
namely, the ratio of incoherent pumping to dissipation in the states
$\kappa/\gamma$ and the ratio $\left({{\omega_{M}}-{\omega_{0}}}\right)/T$. We
analyze how $\kappa$ affects the stationary density matrix of polaritons in
the lower branch (9).

The solution for the stationary density matrix (9) with the coefficients
${P_{N}}$, determined by (12)-(13), allows us to find the average population
of the polaritons $\left\langle{{{\hat{n}}_{j}}}\right\rangle={\mathop{\rm
tr}\nolimits}\left({{{\hat{n}}_{j}}\hat{\rho}}\right)$. At low pumping rates,
the population exponentially decreases with an increase of the state frequency
(Fig. 1). In this case, the exponent does not depend on $\kappa$. However, for
$\kappa$ exceeding a certain threshold value (hereinafter, the condensation
threshold), a relatively large number of polaritons accumulates in the ground
state (Fig. 1). This behavior of the average population of polaritons was
observed in the experiments [6, 10, 19, 13] and associated with the formation
of a polariton BEC. Note that the average polariton density in the ground
state $\left\langle{{{\hat{n}}_{j=0}}}\right\rangle$ sharply increases at the
condensation threshold (Fig. 1) as has been evidenced in various experiments
[6, 10, 19, 13].

Figure 1:  The dependence of the average population of polaritons
$\left\langle{{{\hat{n}}_{j}}}\right\rangle=\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle$
on the incoherent pumping rate $\kappa$ at $M+1=100$,
$\left({{\omega_{M}}-{\omega_{0}}}\right)/T=1$. (Inset) Average population of
polaritons in the ground state.

The incoherent pumping drives the system out of thermodynamic equilibrium.
Therefore, the applicability of equilibrium statistical mechanics approaches
to polariton condensates is at least controversial. To explore differences and
similarities between equilibrium and non-equilibrium BECs, first, we study the
average number of polaritons
$\left\langle{{{\hat{n}}_{j}}}\right\rangle={\rm{Tr}}\left({\hat{\rho}{\hat{n}_{j}}}\right)$
and examine the resultant polariton population with respect to Bose–Einstein
distribution (see Fig. 1)

$\left\langle{{{\hat{n}}_{j}}}\right\rangle=\frac{A}{{\exp\left({{{\left({{\omega_{j}}-{\omega_{0}}-{\mu_{{\rm{eff}}}}}\right)}\mathord{\left/{\vphantom{{\left({{\omega_{j}}-{\omega_{0}}-{\mu_{{\rm{eff}}}}}\right)}{{T_{{\rm{eff}}}}}}}\right.\kern-1.2pt}{{T_{{\rm{eff}}}}}}}\right)-1}}$
(14)

where ${T_{{\rm{eff}}}}$, ${{\mu_{{\rm{eff}}}}}$ and $A$ are extracted by the
least squares fit. For a system in thermodynamic equilibrium
${T_{{\rm{eff}}}}$ is the equilibrium temperature of the system,
${\mu_{{\rm{eff}}}}$ is the chemical potential, and $A$ is an effective
parameter proportional to the number of states around the ground state [54]
defined by size of a system. Like in thermodynamic equilibrium, we assign
${T_{{\rm{eff}}}}$ to the effective condensate temperature,
${{\mu_{{\rm{eff}}}}}$ to the effective chemical potential of the polariton
condensate. Analogously, we suppose that $A$ does not depend on
${T_{{\rm{eff}}}}$ or ${{\mu_{{\rm{eff}}}}}$. Consequently, $A$ should not
depend on the rate of incoherent pumping $\kappa$. The non-equilibrium
polariton condensate demonstrates the same trend of decrease
${\mu_{{\rm{eff}}}}$ with increasing $\kappa$ (Fig. 2), like conventional BECs
in equilibrium statistical mechanics. When the condensation threshold is
exceeded, ${\mu_{{\rm{eff}}}}$ becomes almost equal to zero. However, unlike
BEC at thermodynamic equilibrium the non-equilibrium condensate exhibits
nontrivial dependence of effective temperature ${T_{{\rm{eff}}}}$ on $\kappa$,
as shown in the Figure 2. Indeed, ${T_{{\rm{eff}}}}$ collapses below the
reservoir temperature $T$ above condensation threshold leading to stimulated
cooling of the polariton gas (Fig. 2).

Figure 2:  Approximation of the average population of polaritons in non-
equilibrium BEC by the equilibrium Bose–Einstein distribution (14) at
$M+1=100$, $\left({{\omega_{M}}-{\omega_{0}}}\right)/T=1$. (Blue) The ratio of
the effective temperature of polaritons to the reservoir temperature as a
function of the rate of incoherent pumping. (Red) Effective chemical potential
of polaritons as a function of the rate of incoherent pumping.

Quantum correlations in BEC at the condensation threshold. The main advance
enabled by an exact solution (9) is that it accesses coherence and high-order
correlations of all the modes in polariton BEC. Figure 3 shows the second-
order autocorrelation function for polaritons as a function of $\kappa$. As
$\kappa$ increases, two processes take place: on the one hand, the degree of
coherence increases, and on the other, the frequency region, in which
coherence is formed, narrows. These two processes lead to a non-monotonic
behaviour of the second-order autocorrelation function of the polaritons in
the excited states with respect to the rate of incoherent pumping. At the same
time, the coherence of the ground state monotonically increases. Namely, when
$\kappa$ exceeds the condensation threshold, the second-order coherence
function drops from two to one, see Fig. 3.

Figure 3:  Correlation properties of non-equilibrium polariton BEC. Second-
order autocorrelation function of polaritons in the ground and excited states
versus the rate of incoherent pumping at $M+1=100$,
$\left({{\omega_{M}}-{\omega_{0}}}\right)/T=1$. (Inset) Second-order
autocorrelation function of the ground state.

Another peculiar aspect in the formation of polariton BEC is the correlation
between different polariton states. Although we limit our consideration to
cross-correlations between the ground state and other excited states (Fig. 4),
one can extend the analysis on cross-correlation between arbitrary polariton
states. At the condensation threshold,
$\left\langle{{{\hat{n}}_{0}}{{\hat{n}}_{j\neq
0}}}\right\rangle/\left\langle{{{\hat{n}}_{0}}}\right\rangle\left\langle{{{\hat{n}}_{j\neq
0}}}\right\rangle$ is significantly less than unity. This means that the
ground state and the excited states become correlated in such a way that
occupations of ground and excited states at the same time are incompatible
events.

Figure 4:  Second-order cross-correlation function between the ground
polariton state and excited states versus the rate of incoherent pumping at
$M+1=100$, $\left({{\omega_{M}}-{\omega_{0}}}\right)/T=1$.

Discussion and conclusion In this letter, we report on an analytical quantum
description of non-equilibrium BEC based on solving the master equation for
the full polariton density matrix in the limit of fast polariton
thermalization. While the average polariton density and the effective chemical
potential exhibit the usual properties of BEC at thermal equilibrium, we find
out the effective temperature strongly deviates from the reservoir
temperature, namely polariton gas undergoes noticeable cooling effect. We also
demonstrate that the formation of polariton BEC is accompanied by anti-
correlations between the ground polariton state and excited states gradually
increasing towards the condensation threshold.

It is worth to mention the developed analytical theory strongly relies on the
fast thermalization in bosonic systems. In fact, dynamics of polariton
thermalization depends on the polariton density [14, 23]. Therefore, our
approach is valid for sufficiently intense pumping rate nearby the condensate
threshold [6, 10, 11, 14, 17]. However, for long-lived polariton and photon
systems thermalization occurs at pumping rates significantly lower
condensation threshold [13, 55]. Alongside the conventional strongly-coupled
microcavities we would like to highlight recently developed plasmonic arrays
hybridized with Frenkel excitons as a separate fascinating class of systems
demonstrating BEC phenomena at room temperature with extreme thermalization
rate on the order of 200 fs [56, 57]. The ultra-fast thermalization makes our
theory especially reliable for such systems in general.

We described polaritons as harmonic oscillators (1) thus neglecting polariton
interactions in the system while in some cases, this nonlinearity turns out to
be decisive for the formation of condensate [24, 44]. Overall, the theory
developed in this letter is applicable when the polariton density is large
enough for fast thermalization, but not as high as required for polariton
interactions to show up in the substantial energy change of polariton
eigenstates.

For the sake of simplicity, we considered the case when the incoherent pump
populates the highest-frequency state only. Although in most experiments
pumping configuration is slightly different, there are experimental
realizations analogous to the considered one [8]. Nevertheless, in
Supplemental Material I.6 we provide results obtained under simultaneous
incoherent pump of all considered polariton - the common scenario in
experiments. The dependencies do not exhibit any significant deviation from
the simplified picture considered above.

While the steady-state regime considered in this letter is more illustrative,
and the results are easy to interpret, the general equations (10)–(11) are
applicable to the case of pulsed drive as well. Moreover, recent experimental
advances in quasi-continuous organic polariton condensation [58, 59], bring
this theory into the context of room-temperature polariton BEC.

The work was supported by the Russian Science Foundation (Grant No.
20-72-10145). E.S.A thanks the Foundation for the Advancement of Theoretical
Physics and Mathematics Basis.

## References

  * Bose [1924] S. N. Bose, Zeitschrift für Physik 26, 178 (1924).
  * Einstein [1925] A. Einstein, Sitzungsberichte der Preußischen Akademie der Wissenschaften, Physikalisch-Mathematische Klasse, Berlin , 3 (1925).
  * Anderson  _et al._ [1995] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, science 269, 198 (1995).
  * Davis  _et al._ [1995] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. Kurn, and W. Ketterle, Physical review letters 75, 3969 (1995).
  * Deng  _et al._ [2003] H. Deng, G. Weihs, D. Snoke, J. Bloch, and Y. Yamamoto, Proceedings of the National Academy of Sciences 100, 15318 (2003).
  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. Keeling, F. Marchetti, M. Szymańska, R. André, J. Staehli, _et al._ , Nature 443, 409 (2006).
  * Combescot and Shiau [2015] M. Combescot and S.-Y. Shiau,  _Excitons and Cooper pairs: two composite bosons in many-body physics_ (Oxford University Press, 2015).
  * Byrnes  _et al._ [2014] T. Byrnes, N. Y. Kim, and Y. Yamamoto, Nature Physics 10, 803 (2014).
  * Wertz  _et al._ [2010] E. Wertz, L. Ferrier, D. Solnyshkov, R. Johne, D. Sanvitto, A. Lemaître, I. Sagnes, R. Grousson, A. V. Kavokin, P. Senellart, _et al._ , Nature physics 6, 860 (2010).
  * Balili  _et al._ [2007] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Science 316, 1007 (2007).
  * Estrecho  _et al._ [2018] E. Estrecho, T. Gao, N. Bobrovska, M. D. Fraser, M. Steger, L. Pfeiffer, K. West, T. C. H. Liew, M. Matuszewski, D. W. Snoke, _et al._ , Nature communications 9, 1 (2018).
  * Klaas  _et al._ [2018] M. Klaas, E. Schlottmann, H. Flayac, F. Laussy, F. Gericke, M. Schmidt, M. v. Helversen, J. Beyer, S. Brodbeck, H. Suchomel, _et al._ , Physical review letters 121, 047401 (2018).
  * Sun  _et al._ [2017] Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Physical review letters 118, 016602 (2017).
  * Deng  _et al._ [2010] H. Deng, H. Haug, and Y. Yamamoto, Reviews of Modern Physics 82, 1489 (2010).
  * Imamoglu  _et al._ [1996] A. Imamoglu, R. Ram, S. Pau, Y. Yamamoto, _et al._ , Physical Review A 53, 4250 (1996).
  * Keeling and Kéna-Cohen [2020] J. Keeling and S. Kéna-Cohen, Annual Review of Physical Chemistry 71, 435 (2020).
  * Wei  _et al._ [2019] M. Wei, S. K. Rajendran, H. Ohadi, L. Tropf, M. C. Gather, G. A. Turnbull, and I. D. Samuel, Optica 6, 1124 (2019).
  * Zasedatelev  _et al._ [2019] A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, Nature Photonics 13, 378 (2019).
  * Plumhof  _et al._ [2014] J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt, Nature materials 13, 247 (2014).
  * Kavokin and Malpuech [2003] A. Kavokin and G. Malpuech,  _Thin films and nanostructures_ , Vol. 32 (Academic Press, 2003).
  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen, Nature materials 15, 1061 (2016).
  * Malpuech  _et al._ [2002] G. Malpuech, A. Di Carlo, A. Kavokin, J. J. Baumberg, M. Zamfirescu, and P. Lugli, Applied physics letters 81, 412 (2002).
  * Kavokin  _et al._ [2017] A. Kavokin, J. J. Baumberg, G. Malpuech, and F. P. Laussy,  _Microcavities_ (Oxford university press, 2017).
  * Porras  _et al._ [2002] D. Porras, C. Ciuti, J. Baumberg, and C. Tejedor, Physical Review B 66, 085304 (2002).
  * Banyai and Gartner [2002] L. Banyai and P. Gartner, Physical review letters 88, 210404 (2002).
  * Cao  _et al._ [2004] H. T. Cao, T. Doan, D. T. Thoai, and H. Haug, Physical Review B 69, 245325 (2004).
  * Doan  _et al._ [2008] T. Doan, H. T. Cao, D. T. Thoai, and H. Haug, Physical Review B 78, 205306 (2008).
  * Tassone  _et al._ [1997] F. Tassone, C. Piermarocchi, V. Savona, A. Quattropani, and P. Schwendimann, Physical Review B 56, 7554 (1997).
  * Kirton and Keeling [2013] P. Kirton and J. Keeling, Physical review letters 111, 100404 (2013).
  * Kirton and Keeling [2015] P. Kirton and J. Keeling, Physical Review A 91, 033826 (2015).
  * Strashko  _et al._ [2018] A. Strashko, P. Kirton, and J. Keeling, Physical Review Letters 121, 193601 (2018).
  * Hartwell and Snoke [2010] V. Hartwell and D. Snoke, Physical Review B 82, 075307 (2010).
  * Kasprzak  _et al._ [2008] J. Kasprzak, D. Solnyshkov, R. André, L. S. Dang, and G. Malpuech, Physical review letters 101, 146404 (2008).
  * Sanvitto and Timofeev [2012] D. Sanvitto and V. Timofeev,  _Exciton Polaritons in Microcavities: New Frontiers_ , Vol. 172 (Springer Science & Business Media, 2012).
  * Arnardottir  _et al._ [2020] K. B. Arnardottir, A. J. Moilanen, A. Strashko, P. Törmä, and J. Keeling, Physical Review Letters 125, 233603 (2020).
  * Rubo  _et al._ [2003] Y. G. Rubo, F. Laussy, G. Malpuech, A. Kavokin, and P. Bigenwald, Physical review letters 91, 156403 (2003).
  * Laussy  _et al._ [2004] F. P. Laussy, G. Malpuech, and A. Kavokin, physica status solidi (c) 1, 1339 (2004).
  * del Valle  _et al._ [2009] E. del Valle, D. Sanvitto, A. Amo, F. P. Laussy, R. André, C. Tejedor, and L. Viña, Phys. Rev. Lett. 103, 096404 (2009).
  * Rubo [2004] Y. G. Rubo, physica status solidi (a) 201, 641 (2004).
  * Laussy  _et al._ [2003] F. P. Laussy, Y. G. Rubo, G. Malpuech, A. Kavokin, and P. Bigenwald, physica status solidi (c) n/a, 1476 (2003).
  * Lagoudakis  _et al._ [2003] P. Lagoudakis, M. Martin, J. Baumberg, A. Qarry, E. Cohen, and L. Pfeiffer, Physical review letters 90, 206401 (2003).
  * Maragkou  _et al._ [2010] M. Maragkou, A. Grundy, T. Ostatnickỳ, and P. Lagoudakis, Applied Physics Letters 97, 111110 (2010).
  * Coles  _et al._ [2011] D. M. Coles, P. Michetti, C. Clark, W. C. Tsoi, A. M. Adawi, J.-S. Kim, and D. G. Lidzey, Advanced Functional Materials 21, 3691 (2011).
  * Savvidis  _et al._ [2000] P. Savvidis, J. Baumberg, R. Stevenson, M. Skolnick, D. Whittaker, and J. Roberts, Physical review letters 84, 1547 (2000).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Scully and Zubairy [1999] M. O. Scully and M. S. Zubairy, Quantum optics (1999).
  * Litinskaya  _et al._ [2004] M. Litinskaya, P. Reineker, and V. Agranovich, Journal of luminescence 110, 364 (2004).
  * Mazza  _et al._ [2009] L. Mazza, L. Fontanesi, and G. La Rocca, Physical Review B 80, 235314 (2009).
  * Bittner and Silva [2012] E. R. Bittner and C. Silva, The Journal of chemical physics 136, 034510 (2012).
  * Ćwik  _et al._ [2014] J. A. Ćwik, S. Reja, P. B. Littlewood, and J. Keeling, EPL (Europhysics Letters) 105, 47009 (2014).
  * Somaschi  _et al._ [2011] N. Somaschi, L. Mouchliadis, D. Coles, I. Perakis, D. Lidzey, P. Lagoudakis, and P. Savvidis, Applied Physics Letters 99, 209 (2011).
  * Ramezani  _et al._ [2018] M. Ramezani, Q. Le-Van, A. Halpin, and J. G. Rivas, Physical Review Letters 121, 243904 (2018).
  * Shishkov  _et al._ [2018] V. Y. Shishkov, E. Andrianov, A. Pukhov, A. Vinogradov, and A. Lisyansky, Physical Review E 98, 022132 (2018).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Statistical Physics: Volume 5_ , Vol. 5 (Elsevier, 2013).
  * Weill  _et al._ [2019] R. Weill, A. Bekker, B. Levit, and B. Fischer, Nature communications 10, 1 (2019).
  * Hakala  _et al._ [2018] T. K. Hakala, A. J. Moilanen, A. I. Väkeväinen, R. Guo, J.-P. Martikainen, K. S. Daskalakis, H. T. Rekola, A. Julku, and P. Törmä, Nature Physics 14, 739 (2018).
  * Väkeväinen  _et al._ [2020] A. I. Väkeväinen, A. J. Moilanen, M. Nečada, T. K. Hakala, K. S. Daskalakis, and P. Törmä, Nature communications 11, 1 (2020).
  * Betzold  _et al._ [2019] S. Betzold, M. Dusel, O. Kyriienko, C. P. Dietrich, S. Klembt, J. Ohmer, U. Fischer, I. A. Shelykh, C. Schneider, and S. Höfling, ACS Photonics 7, 384 (2019).
  * Putintsev  _et al._ [2020] A. Putintsev, A. Zasedatelev, K. E. McGhee, T. Cookson, K. Georgiou, D. Sannikov, D. G. Lidzey, and P. Lagoudakis, Applied Physics Letters 117, 123302 (2020).

##  I SUPPLEMENTAL MATERIAL

###  I.1  Properties of the partition function ${Z_{N}}$

Some important properties of the partition function ${Z_{N}}$, defined by (8),
can be obtained using its generating function

$f\left(x\right)=\sum\limits_{N=0}^{+\infty}{{x^{N}}{Z_{N}}}=\prod\limits_{j=0}^{M}{\frac{1}{1-{w_{j}}x}},$
(15)

where $M+1$ is the total number of polariton states,
${w_{j}}=\exp\left({-\hbar\left({{\omega_{j}}-{\omega_{0}}}\right)/T}\right)$,
${\omega_{j}}$ is the frequency of the $j$th state and $T$ is the reservoir
temperature. The generating function $f\left(x\right)$ is defined for $0\leq
x<1$.

Using the generating function (15), one can obtain the following properties of
${Z_{N}}$

$\mathop{\lim}\limits_{N\to+\infty}{Z_{N}}=\mathop{\lim}\limits_{x\to
1-0}\left({1-x}\right)f\left(x\right)$ (16)
$w_{j}\frac{\partial{Z_{N}}}{\partial{w_{j}}}=\sum\limits_{n=0}^{N}{w_{j}^{n}{Z_{N-n}}}-Z_{N}$
(17)

###  I.2 Derivation of Eqs. (10)-(11)

Substitution of the general expression for the density matrix in the limit of
fast thermalization (9) into the master equation (5) leads to

$\sum\limits_{N=0}^{+\infty}{\frac{\partial{P_{N}}\left(t\right)}{\partial
t}}{\frac{1}{{Z_{N}}}}\sum\limits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}{{\hat{r}}_{{n_{0}}...{n_{M}}}}}=\\\
=\sum\limits_{N=0}^{+\infty}\frac{{P_{N}}\left(t\right)}{Z_{N}}\sum\limits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}{{L}_{{\rm{diss}}}}\left({{{\hat{r}}_{{n_{0}}...{n_{M}}}}}\right)}+\sum\limits_{N=0}^{+\infty}\frac{{P_{N}}\left(t\right)}{Z_{N}}\sum\limits_{{n_{0}}+...+{n_{M}}=N}{w_{0}^{{n_{0}}}...w_{M}^{{n_{M}}}{{L}_{{\rm{pump}}}}\left({{{\hat{r}}_{{n_{0}}...{n_{M}}}}}\right)}$
(18)

where
${\hat{r}_{{n_{0}}...{n_{M}}}}=\left|{{n_{0}},...,{n_{M}}}\right\rangle\left\langle{{n_{0}},...,{n_{M}}}\right|$.
To obtain Equations (10)-(11) it is necessary to project Equation (18) to
states with the total number of polaritons $N$.

###  I.3 Density matrix of the polaritons in the ground state

Figure 5 shows the diagonal elements of the polariton density matrix
$\rho_{N,N}^{{\rm{condensate}}}$ for $j=0$, defined by
${\hat{\rho}_{0}}=\sum\nolimits_{N}{\left|N\right\rangle\left\langle
N\right|\rho_{N,N}^{{\rm{condensate}}}}$, where
${\hat{\rho}_{0}}={{\mathop{\rm tr}\nolimits}_{1}}{{\mathop{\rm
tr}\nolimits}_{2}}...{{\mathop{\rm
tr}\nolimits}_{M}}\left({\hat{\rho}}\right)$, $\hat{\rho}$ is the full density
matrix of the lower polaritons and ${{\mathop{\rm tr}\nolimits}_{j}}$ is
partial trace over the $j$th state. As one can see from Figure 5, at pumping
rates below the threshold, the maximal diagonal element of the density matrix
is $\rho_{0,0}^{{\rm{condensate}}}$. However, at pumping rates above the
threshold, the maximal diagonal element of the density matrix
$\rho_{N,N}^{{\rm{condensate}}}$ corresponds to $N>0$, and the shape of the
distribution of the diagonal elements of the density matrix is getting close
to a coherent distribution.

Figure 5:  Formation of BEC. Diagonal elements of the density matrix of the
ground state at different rates of incoherent pumping.

###  I.4 Approximation of the average population of polaritons by the
Bose–Einstein distribution

Figure 6 shows the calculated distribution of the average population of
polaritons and its approximation by the Bose–Einstein distribution. Figure 6
shows the relative error arising from this approximation. One can see that the
Bose–Einstein distribution approximates quite well the stationary distribution
of the non-equilibrium polariton BEC.

Figure 6:  (a) Calculated average population of polaritons. (b) Approximation
of the average population of polaritons by the Bose–Einstein distribution with
parameter $A=1.9$ (see (14)). (c) The relative error between the calculated
population and the Bose–Einstein distribution. Figure 7:  A decrease in the
effective temperature of polaritons with incoherent pumping rates above the
threshold of the BEC formation.

###  I.5 Dynamics of the formation of non-equilibrium BEC.

Figure 8 shows the time dependence of the population of polaritons in the
ground state above the threshold. One can see that the build-up of a
stationary population of the ground state occurs in a characteristic time
equal to the decay time of polaritons. In addition, the higher the rate of the
incoherent pumping, the faster the build-up of the stationary population in
the ground state.

Figure 8:  The ratio of the steady-state population of the ground state and
current population of the ground state for pumping rates above the condensate
formation threshold.

###  I.6 Formation of polariton BEC when all states of the lower polariton
branch are pumped

Figure 9 shows the dependence of the average values of polariton population at
different pumping rates. In this case, incoherent pumping acts on all the
polariton states in such a way that its effect on the polaritons can be
described by the Lindblad operator (3) from the main text with
${\kappa_{j}}=\kappa$ for all polariton states. It can be seen from the figure
that there are no qualitative differences for the stationary average
population of polaritons between the cases when only the highest-frequency
polariton state is pumped and when all polariton states are pumped
simultaneously.

Figure 9:  The average population of lower polaritons when all states are
pumped simultaneously.