#  Analytical framework for non-equilibrium phase transition to Bose–Einstein
condensate

V. Yu. Shishkov<EMAIL_ADDRESS>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;  Center for Photonics and Quantum Materials, Skolkovo
Institute of Science and Technology, Moscow, Russian Federation  Laboratories
for Hybrid Photonics, Skolkovo Institute of Science and Technology, Moscow,
Russian Federation  E. 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;  Center for Photonics and Quantum Materials, Skolkovo Institute of
Science and Technology, Moscow, Russian Federation  Laboratories for Hybrid
Photonics, Skolkovo Institute of Science and Technology, Moscow, Russian
Federation  Yu. E. Lozovik  Dukhov Research Institute of Automatics (VNIIA),
22 Sushchevskaya, Moscow 127055, Russia;  Moscow Institute of Electronics and
Mathematics, National Research University Higher School of Economics, 101000
Moscow, Russia;  Institute for Spectroscopy RAS, 5 Fizicheskaya, Troitsk
142190, Russia;

###### Abstract

The theoretical description of non-equilibrium Bose–Einstein condensate (BEC)
is one of the main challenges in modern statistical physics and kinetics. The
non-equilibrium nature of BEC makes it impossible to employ the well-
established formalism of statistical mechanics. We develop a framework for the
analytical description of a non-equilibrium phase transition to BEC that, in
contrast to previously developed approaches, takes into account the infinite
number of continuously distributed states. We consider the limit of fast
thermalization and obtain an analytical expression for the full density matrix
of a non-equilibrium ideal BEC which also covers the equilibrium case. For the
particular cases of 2D and 3D, we investigate the non-equilibrium formation of
BEC by finding the temperature dependence of the ground state occupation and
second-order coherence function. We show that for a given pumping rate, the
macroscopic occupation of the ground state and buildup of coherence may occur
at different temperatures. Moreover, the buildup of coherence strongly depends
on the pumping scheme. We also investigate the condensate linewidth and show
that the Schawlow–Townes law holds for BEC in 3D and does not hold for BEC in
2D.

##  1 Introduction

Non-equilibrium Bose–Einstein condensates (BECs) working at room temperature
[1, 2] are systems where open dissipative quantum-mechanical nature meets
many-body collective dynamics. Besides its fundamental theoretical interest,
BEC can be used in sub-picosecond switching at the fundamental quantum limit
[3], all-optical manipulation, and other low-energy optoelectronic
applications [4, 5]. Among the most promising systems for the realization of
BEC are exciton-polaritons, the hybrid quasiparticles formed by the strong
light–matter interaction of photons in a cavity and excitons, [6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 5, 17, 2, 1]. The advantage of the exciton-polaritons
is their small effective mass and low density of states [9, 10, 14, 16], that
provides a critical temperature of hundreds of kelvin [6, 7, 1].

In many BEC realizations, the thermalization rate of the polaritons is much
faster than their dissipation rate and pumping rate. Namely, for long-lived
polariton and photon systems, thermalization occurs at pumping rates
significantly lower than the condensation threshold [13, 18]. The relatively
high thermal energy in polariton systems at room temperature facilitates
phonon- and vibron-assisted polariton relaxation, which, overall, results in
efficient polariton thermalization. A fast thermalization in an organic
polariton system has been revealed with an extreme rate on the order of 200 fs
[19, 20].

Previous studies of BEC in these systems largely rely on the semiclassical
Maxwell–Boltzmann equations describing the average population of exciton-
polariton states [21, 22, 23, 24, 25], other mean-field theories [26, 27, 28],
and an approach beyond mean-field theory [29] that takes into account the
higher-order correlations. However, in such treatments, information about the
correlations between different polariton states is lost. In many cases, the
driven dissipative Gross–Pitaevskii equation with noise successfully describes
the dynamics of BEC [30]. However, the Gross–Pitaevskii equation works well
only when temperature is well below the condensation temperature [30]. This
means that it cannot accurately describe the crossover and buildup of
coherence.

The dynamics of the polaritons can be described by the corresponding master
equation for the density matrix [31, 32], with transition rates between
quantum states in a configuration space rather than between populations as in
semiclassical Maxwell–Boltzmann equations. The density matrix enables finding
correlations of any order, but the master equation become hard to solve
numerically for large systems. The analytical steady state solution for the
density matrix for two-oscillators system were obtained in [33, 34] and
generalized for arbitrary degeneracy of the excited state in [32]. In these
works, the authors discussed how the number-conserving scattering processes
allow to obtain this analytical solution. Relying on similar ideas, recently,
it has been shown that, in the case of fast thermalization of polaritons, the
density matrix can be obtained analytically for arbitrary finite sets of
polariton states [35]. However, due to the finiteness of the states, this
approach is limited in its ability to analytically investigate the influence
of the polaritons’ dispersion and dimension on the condensation threshold,
buildup of coherence, long-rage spatial correlations, and linewidth. Thus, a
comprehensive analysis of BEC requires the extension of the approach [35] to
an infinite number of states.

In the present paper, we develop a formalism which extends the advantages of
the exactly solvable Lindblad master equation for non-equilibrium BEC in the
fast thermalization limit [35] to the systems with continuously distributed
states. This developed approach consists of two essential steps. At the first
step, it is assumed that in each sector of Fock space with fixed particle
number $N$ distributed in the whole system, the density matrix
$\hat{\rho}_{N}$ is in thermal equilibrium. This assumption is fulfilled if
thermalizing processes are much faster than processes that change the particle
number (external pumping and dissipation) [35]. In this case, the total system
density matrix can be written as a time-dependent sum of thermal, canonical
density matrices in each particle number sector. Thus, the problem is
effectively reduced to calculating all thermal density matrices for arbitrary,
but fixed particle number. At the second step, we develop a recursive method
for calculating all equilibrium density matrices with fixed particle number.
As a result, we obtain the full density matrix of all the states. We show
that, for a non-equilibrium BEC in the fast thermalization limit, the
condition for the macroscopic occupation of the ground state is the same as in
the full thermal equilibrium case. However, the buildup of coherence in a non-
equilibrium BEC depends not only on the number of polaritons and the
environmental temperature, but also on the pumping scheme. An increase in the
pumping rate of the ground state lead to the macroscopic occupation of the
ground state, but destroys the coherence of a non-equilibrium BEC; on the
contrary, an increase in the pumping rate of the excited states results in
both the macroscopic occupation of the ground state and the buildup of
coherence. We also show that the linewidth of the condensate narrows with an
increase of the macroscopic occupation of the ground state, but in 2D this
narrowing does not follow the Schawlow–Townes equation.

##  2 Description of the model

The strong light–matter interaction of a cavity electric field with the dipole
moments of an optical transition of molecules gives rise to new eigenstates,
namely, lower and upper polariton branches. Due to dissipation processes, it
is necessary to compensate for the losses, which can be done by an external
pump. Usually, an external pump excites polaritons in the upper polariton
branch and excitons which are not coupled with the cavity modes [2, 9]. The
polaritons of the upper branch and excitons [2, 9] scatter the energy supplied
by the pump into the lower polariton branch forming BEC at high enough
excitation density. There are several mechanisms allowing such energy
relaxation in polariton systems. Scattering by electrons [36], high-energy
optical phonons [37], and vibrons in organic materials [2, 38] are among them.
Since non-equilibrium BEC occurs mostly 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
by an effective incoherent pumping. Such an approximation is reasonable for
the majority of experimental realizations [7, 1, 13]

We consider a continuous system of polariton states with the frequency
$\omega_{\bf k}$ for a wave vector ${\bf k}$ in the general case of $D$
dimensional space. We suppose that the polaritons with wave vector ${\bf k}$
of the lower polariton branch are described by bosonic creation $\hat{a}_{\bf
k}^{\dagger}$ and annihilation $\hat{a}_{\bf k}$ operators [31, 32] which obey
the commutation relation $\left[{{{\hat{a}}_{\bf k}},\hat{a}_{\bf
q}^{\dagger}}\right]={\delta({\bf k}-{\bf q})}$. In this case, the Hamiltonian
of the polaritons takes the form

${\hat{H}_{{\rm{LP}}}}=\sum_{\bf k}{{\omega_{\bf k}}\hat{a}_{\bf
k}^{\dagger}{{\hat{a}}_{\bf k}}},$ (1)

where the state with ${\bf k}={\bf 0}$ is the ground state. Hereafter, we
assume that $\hbar=1$. We consider the ideal polaritons, i.e. there is no
interaction between polaritons in the Hermitian part.

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

${L_{{\rm{diss}}}}\left({\hat{\rho}}\right)=\sum_{\bf k}{{\gamma_{\bf
k}}\left({{{\hat{a}}_{\bf k}}\hat{\rho}\hat{a}_{\bf
k}^{\dagger}-\frac{1}{2}\hat{\rho}\hat{a}_{\bf k}^{\dagger}{{\hat{a}}_{\bf
k}}-\frac{1}{2}\hat{a}_{\bf k}^{\dagger}{{\hat{a}}_{\bf
k}}\hat{\rho}}\right)},$ (2)

where ${\gamma_{\bf k}}$ is the dissipation rate of the state with wave vector
${\bf k}$. The effective incoherent pumping of the lower polaritons can be
described by the Lindblad superoperator

${L_{{\rm{pump}}}}\left({\hat{\rho}}\right)=\\\ \sum_{\bf k}{{\kappa_{\bf
k}}\left({{{\hat{a}}_{\bf k}}\hat{\rho}\hat{a}_{\bf
k}^{\dagger}-\frac{1}{2}\hat{\rho}\hat{a}_{\bf k}^{\dagger}{{\hat{a}}_{\bf
k}}-\frac{1}{2}\hat{a}_{\bf k}^{\dagger}{{\hat{a}}_{\bf
k}}\hat{\rho}}\right)}\\\ +\sum_{\bf k}{{\kappa_{\bf k}}\left({\hat{a}_{\bf
k}^{\dagger}\hat{\rho}{{\hat{a}}_{\bf k}}-\frac{1}{2}\hat{\rho}{{\hat{a}}_{\bf
k}}\hat{a}_{\bf k}^{\dagger}-\frac{1}{2}{{\hat{a}}_{\bf k}}\hat{a}_{\bf
k}^{\dagger}\hat{\rho}}\right)},$ (3)

where ${\kappa_{\bf k}}$ is the pumping rate of the state with wave vector
${\bf k}$. The Lindblad operator $L_{\rm pump}$ leads to the following
dynamics of the average polariton population in the state with the wave vector
$\bf k$ [35]: $\left(d\langle\hat{a}_{\bf k}^{\dagger}\hat{a}_{\bf
k}\rangle/dt\right)_{\rm pump}=\kappa_{\bf k}$ and $\left(d\langle\hat{a}_{\bf
k}\rangle/dt\right)_{\rm pump}=0$. Therefore, $L_{\rm pump}$ leads to the
excitation of a certain number of polaritons per unit time in the
corresponding state without affecting the phase of the state.

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 [28, 40, 41, 42, 43, 44, 45]. For polariton states in inorganic
semiconductors, thermalization predominantly goes through interactions with
acoustic phonons or free charges [31]. The thermalization of the polaritons
can also occur due to polariton–polariton scattering [14, 46]. However,
regardless of the mechanism, the thermalization can be described through the
Lindblad superoperator [31, 32]

${L_{{\rm{therm}}}}\left({\hat{\rho}}\right)=\sum_{\bf k}\sum_{{\bf q}\neq{\bf
k}}{\Gamma_{\bf kq}}\times\\\ {\left({{{\hat{a}}_{\bf k}}\hat{a}_{\bf
q}^{\dagger}\hat{\rho}{{\hat{a}}_{\bf q}}\hat{a}_{\bf
k}^{\dagger}-\frac{1}{2}\hat{\rho}{{\hat{a}}_{\bf q}}\hat{a}_{\bf
k}^{\dagger}{{\hat{a}}_{\bf k}}\hat{a}_{\bf
q}^{\dagger}-\frac{1}{2}{{\hat{a}}_{\bf q}}\hat{a}_{\bf
k}^{\dagger}{{\hat{a}}_{\bf k}}\hat{a}_{\bf q}^{\dagger}\hat{\rho}}\right)},$
(4)

where ${\Gamma_{\bf kq}}$ is the transition rate from the state with the wave
vector $\bf k$ to the state with the wave vector $\bf q$. The thermalization
rates ${\Gamma_{\bf kq}}$ obey the Kubo–Martin–Schwinger relation
${\Gamma_{\bf kq}}/{\Gamma_{\bf qk}}=\exp\left({\left({{\omega_{\bf
k}}-{\omega_{\bf q}}}\right)/T}\right)$, where $T$ is the temperature of
intermolecular vibrations of the organic dyes or the temperature of the
phonons in the semiconductors [47]. Hereafter, we assume the Boltzmann
constant to be equal to unity.

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

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

##  3 Density matrix in the fast thermalization limit

Further, we assume that the thermalization is the fastest process in the
system, i.e., thermalization rates are much higher than the dissipation and
pumping rates, ${\Gamma_{\bf 0k}}(1+\langle\hat{a}^{\dagger}_{\bf
0}\hat{a}_{\bf 0}\rangle)\gg{\kappa_{\bf k}},{\gamma_{\bf k}}$. In such a
case, at times $t\ll\gamma_{\bf k}^{-1},\kappa_{\bf k}^{-1}$, only
thermalization affects the system dynamics. As a result, the density matrix
obeys the approximate differential equation

$\frac{{\partial\hat{\rho}(t)}}{{\partial
t}}=-i\left[{\hat{\rho}(t),{{\hat{H}}_{{\rm{LP}}}}}\right]+{L_{{\rm{therm}}}}\left({\hat{\rho}(t)}\right).$
(6)

After the time ${\Gamma_{\bf 0k}}^{-1}(1+\langle\hat{a}^{\dagger}_{\bf
0}\hat{a}_{\bf 0}\rangle)^{-1}\ll t\ll\gamma_{\bf k}^{-1},\kappa_{\bf
k}^{-1}$, the system reaches its quasistationary state determined by

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

The relaxation operator ${L_{{\rm{therm}}}}\left({\hat{\rho}}\right)$, defined
by the expression (4), conserves the total number of lower polaritons [31].
Thus, in the limit ${\kappa_{\bf k}}/{\Gamma_{\bf
0k}}(1+\langle\hat{a}^{\dagger}_{\bf 0}\hat{a}_{\bf 0}\rangle)\to 0$ and
${\gamma_{\bf k}}/{\Gamma_{\bf 0k}}(1+\langle\hat{a}^{\dagger}_{\bf
0}\hat{a}_{\bf 0}\rangle)\to 0$, the system reaches a quasistationary state
governed by Eq. (7) at instantaneous time conserving the total number of
polaritons. Applying the theory developed in [48], we can obtain the general
form of the density matrix of this quasistationary state. Because the
relaxation operator ${L_{{\rm{therm}}}}\left({\hat{\rho}}\right)$ conserves
the total number of lower polaritons [31], the operator of total number of
lower polaritons $\sum_{\bf k}\hat{a}^{\dagger}_{\bf k}\hat{a}_{\bf k}$ is an
integral of motion for the thermalization process. According to [48], the
presence of this integral of motion implies that the system has invariant
subspaces with the total number of polaritons equal to $\sum_{\bf k}{n_{\bf
k}}=N$. Being in an invariant subspace at the initial moment in time, the
system stays in this invariant subspace in the subsequent moments. In each
invariant subspace with $N$ polaritons, the Gibbs distribution with a
temperature $T$ is established [48].

The most general form of density matrix satisfying Eq. (7) is the weighted sum
of thermalized distributions in each of the invariant subspaces:

$\hat{\rho}(t)=\sum_{N=0}^{+\infty}P_{N}(t)\hat{\rho}_{N},$ (8)

where the summation goes over the total number of polaritons. Here $P_{N}(t)$
is the probability of the system to be in a state with exactly $N$ polaritons
distributed over the system as a whole. These probabilities depend on the time
due to the dissipation and pumping processes. ${\hat{\rho}}_{N}$ is the
thermalized density matrix of the states with exactly $N$ polaritons, which
forms an invariant subspace [35]. This means that $\hat{\rho}_{N}$ can be
represented as

${\hat{\rho}_{N}}=\frac{1}{{{Z_{N}}}}\sum\limits_{{\rm
config}:\sum\limits_{\bf{k}}{{n_{\bf{k}}}}=N}{{e^{-\sum\limits_{\bf{k}}{{n_{\bf{k}}}}({\omega_{\bf{k}}}-{\omega_{\bf{0}}})/T}}}{\hat{R}_{\rm{config}}},$
(9)

where $Z_{N}$ is the partition function of the states with exactly $N$
polaritons, and ${\hat{R}_{\rm{config}}}$ is a diagonal density matrix
corresponding to the configuration of the polaritons $\\{n_{\bf k}\\}$. This
density matrix obeys ${{\hat{a}^{\dagger}_{\bf k}}{\hat{a}_{\bf
k}}}{\hat{R}_{{\rm config}}}=n_{\bf k}{\hat{R}}_{\rm{config}}$,
$\sum_{\bf{k}}{{n_{\bf{k}}}}=N$. We note, that according to [48] the matrix
(9) is diagonal. The consideration of diagonal elements is enough to
investigate many important properties of BEC such as stationary density
matrix, occupation numbers, entropy, coherence, however, we do take into
account the dynamics of non-diagonal elements while finding spectral linewidth
in Section 5.

The partition function $Z_{N}$ is defined by the normalization condition
$\rm{Tr}({\hat{R}}_{\rm{config}})=1$, thus

$Z_{N}=\sum_{{\rm config}:\sum_{\bf k}n_{\bf k}=N}e^{-\sum_{\bf k}n_{\bf
k}(\omega_{\bf k}-\omega_{\bf 0})/{T}}.$ (10)

When there are no polaritons in the system, Eq. (9) implies
${\hat{\rho}}_{N=0}=Z_{0}^{-1}|0\rangle\langle 0|$. Therefore,

$Z_{0}=1.$ (11)

The density matrix (8) allows to calculate all the single-time averages. In
particular, it allows calculating the averages $\langle\hat{n}_{{\bf
k}_{1}}^{m_{1}}\leavevmode\nobreak\ ...\leavevmode\nobreak\ \hat{n}_{{\bf
k}_{M}}^{m_{M}}\rangle$:

$\langle\hat{n}_{{\bf k}_{1}}^{m_{1}}\leavevmode\nobreak\
...\leavevmode\nobreak\ \hat{n}_{{\bf k}_{M}}^{m_{M}}\rangle=\\\
\left(-T\right)^{\sum_{j=1}^{M}m_{j}}\sum_{N=0}^{+\infty}\frac{P_{N}(t)}{Z_{N}}\frac{{\partial^{m_{1}}}\leavevmode\nobreak\
...\leavevmode\nobreak\ {\partial^{m_{M}}}Z_{N}}{{\partial\omega_{{\bf
k}_{1}}^{m_{1}}}\leavevmode\nobreak\ ...\leavevmode\nobreak\
{\partial\omega_{{\bf k}_{M}}^{m_{M}}}}.$ (12)

The derivative of $Z_{N}$ with respect to $\omega_{\bf k}$ is (see Appendix A)

$\frac{\partial Z_{N}}{\partial\omega_{\bf
k}}=\frac{1}{T}\left[Z_{N}-\sum_{n=0}^{N}Z_{N-n}e^{-n(\omega_{\bf
k}-\omega_{\bf 0})/{T}}\right].$ (13)

Note that $\partial Z_{N}/\partial\omega_{\bf 0}$ must be calculated as
$\lim_{{\bf k}\to{\bf 0}}\partial Z_{N}/\partial\omega_{\bf k}$.

##  4 Partition function $Z_{N}$

The aim of paper is to study the process of formation of BEC and its
properties in the continuous limit. This limit covers not only the system
occupying the infinite volume, but also the finite volume as well. Indeed, the
states may be treated as a continuum for systems of finite size when the
volume is so large that the energy spacing between the quantized states is
small compared to $T$ [49]. However, even if this condition is not satisfied,
the analysis presented below may be useful for qualitative analysis of BEC.

We consider the continuous limit of the partition function (10) and replace
the sum over wave vectors $\bf k$ by the integral over frequencies $\omega$:

$\sum_{\bf k}\leavevmode\nobreak\ ...=\int
d\omega\left(1\cdot\delta(\omega-\omega_{\bf
0})+\nu(\omega)\right)\leavevmode\nobreak\ ...$ (14)

The density of states $\nu(\omega)$ is

$\nu(\omega)=\frac{V}{(2\pi)^{D}}\left(\frac{d^{D}{\bf k}}{d\omega_{\bf
k}}\right)_{\omega_{\bf k}=\omega},$ (15)

where $V$ is the volume in 3D or the area in 2D occupied by polaritons. The
delta function with factor unity in the right-hand side of Eq. (14) suggests
that ${\bf k}={\bf 0}$ is non-degenerate. The separation of the state ${\bf
k}={\bf 0}$ is a standard step in the description of continuous systems [50,
51, 31]. Below we consider the quadratic dispersion, $\omega_{\bf
k}=\omega_{\bf 0}+\alpha{\bf k}^{2}$. This dispersion is characteristic for
polaritons of the lower branch in the vicinity of ${\bf k}={\bf 0}$ and the
most frequent case in the experimental realizations of BEC.

The definition of the partition function (10) cannot be used directly to
calculate $Z_{N}$ in the continuous limit. The fastest way to calculate
$Z_{N}$ is to use a recurrence relation [52, 53, 54, 55, 56, 57]. The
recurrence relation follows from the fact that the density matrix (9)
corresponds to exactly $N$ polaritons distributed in the system as a whole.
Thus,

$N=\sum_{\bf k}\langle\hat{n}_{\bf k}\rangle_{N},$ (16)

where $\langle\hat{n}_{\bf k}\rangle_{N}$ stands for ${\rm Tr}(\hat{n}_{\bf
k}\hat{\rho}_{N})$. According to Eqs. (12) and (13),

$\sum_{\bf k}\langle\hat{n}_{\bf
k}\rangle_{N}=\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}\sum_{\bf k}e^{-n(\omega_{\bf
k}-\omega_{\bf 0})/{T}}.$ (17)

After integration over $\bf k$ we obtain the recurrence relations

$Z_{N}=\frac{1}{N}\sum_{n=0}^{N-1}Z_{n}\left(1+\frac{G_{2D}}{N-n}\right)$ (18)

with $G_{2D}=VT/4\pi\alpha$ in 2D and

$Z_{N}=\frac{1}{N}\sum_{n=0}^{N-1}Z_{n}\left(1+\frac{G_{3D}}{(N-n)^{3/2}}\right)$
(19)

with $G_{3D}=V(T/4\pi\alpha)^{3/2}$ in 3D. In both cases, the initial value is
$Z_{0}=1$ (see Eq. (11)). The physical meaning of $G_{2D}$ and $G_{3D}$ is the
number of states in the energy range $T$ above ${\bf k}={\bf 0}$ in the
corresponding dimensions.

The exact expressions for $Z_{N}$ are obtained in Appendix A trough
Darwin–Fowler method [58, 52, 59, 60, 61]. However, for practical
calculations, the optimal strategy to calculate $Z_{N}$ is to use the
recurrence relations (18) and (19).

We can distinguish three regions in $\\{N,T\\}$ plane and write the
asymptotics for $Z_{N}$ in the corresponding regimes.

#### BEC is not formed.

This regime is defined by $1\ll N\ll G_{2D}$ in 2D and $1\ll N\ll G_{3D}$ in
3D ($1\ll N\ll T/T_{0}$ in the figures). From the Eqs. (18) and (19) we obtain

$Z_{N}\approx\frac{G_{2D}^{N}}{N!}$ (20)

in 2D and

$Z_{N}\approx\frac{G_{3D}^{N}}{N!}$ (21)

in 3D. These asymptotics corresponds to the ideal Boltzmann gas. This is
expected, since in this region the concentration of polaritons is small,
therefore, the Bose effect does not have an influence on the distribution of
the polaritons.

#### BEC is formed.

This regime is defined by $1\ll G_{2D}\ll N$ in 2D and $1\ll G_{3D}\ll N$ in
3D ($1\ll T/T_{0}\ll N$ in the figures). From the Eqs. (80) and (82) we obtain

$\sum_{n=0}^{N-1}Z_{n}\approx Z_{\infty}(N-G_{2D}\ln N)\;\Rightarrow\;\\\
\;\Rightarrow\;Z_{N}\approx Z_{\infty}\left(1-\frac{G_{2D}}{N}\right)$ (22)

with $Z_{\infty}=\exp(G_{2D}\zeta(2))$ in 2D and

$\sum_{n=0}^{N-1}Z_{n}\approx Z_{\infty}(N-G_{3D}\zeta(3/2))\;\Rightarrow\;\\\
\;\Rightarrow\;Z_{N}\approx Z_{\infty}\left(1-\frac{G_{3D}}{N^{3/2}}\right)$
(23)

with $Z_{\infty}=\exp(G_{3D}\zeta(5/2))$ in 3D, where $\zeta(x)$ is the zeta
function evaluated at $x$. This regime corresponds to the formation of BEC.

#### Deep quantum regime.

This regime is defined by $G_{2D}\ll 1$ in 2D and $G_{3D}\ll 1$ ($T/T_{0}\ll
1$ in the figures), when almost only ground state is populated. From the Eqs.
(18) and (19) both in 2D and in 3D we obtain

$Z_{N}\approx 1.$ (24)

##  5 Condensation of exactly $N$ polaritons

Figure 1:  The fraction of the polaritons in the excited states in 2D (a) and
in 3D (c), given that there are exactly $N$ polaritons distributed in the
system as a whole. The entropy per polariton in 2D (b) and in 3D (d), given
that there are exactly $N$ polaritons in the system. $T$ is the temperature of
the environment, $T_{0}$ is the temperature corresponding to $G_{2D}=1$ in 2D
and $G_{3D}=1$ in 3D (see Section 4). The black dashed line marks the
approximate condition for a macroscopic occupation of the ground state
$N=G_{2D}\ln N$ in 2D (see Eq. (26)) and $N=G_{3D}\zeta(3/2)$ in 3D (see Eq.
(27)).

In Section 3 we showed that in the fast thermalization limit the invariant
subspaces play a substantial role in the dynamics of BEC. In this section, we
explore the formation of BEC in these invariant subspaces. This means that in
Eq. (8) only one probability $P_{N}$ is non-zero. Therefore, the system is
described by the density matrix (9) and has exactly $N$ polaritons distributed
in the system as a whole.

###  5.1 Population of the states

From Eq. (12) it follows that the average number of polaritons in the state
with the wave vector $\bf k$ can be calculated as

$\langle\hat{n}_{\bf
k}\rangle_{N}=\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}e^{-n(\omega_{\bf
k}-\omega_{\bf 0})/{T}}.$ (25)

The subscript $N$ refers to the total number of polaritons distributed in the
system as a whole. This expression is the same both in 2D and in 3D. From the
Fig. 1 one can see that the macroscopic occupation occurs for the fixed volume
both in 2D and in 3D at sufficiently low temperatures.

In the limit $N\to\infty$, from Eq. (25) we obtain the fraction of polaritons
in the ground state (see Appendix B)

$\frac{\langle\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}}{N}\approx
1-\left(\frac{T}{4\pi\alpha}\right)\left(\frac{V}{N}\right)\ln(N)$ (26)

in 2D and

$\frac{\langle\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}}{N}\approx
1-\left(\frac{T}{4\pi\alpha}\right)^{3/2}\left(\frac{V}{N}\right)\zeta(3/2)$
(27)

in 3D.

In the limit $N,V\to\infty$, $N/V\to\rm{const}$, the Eq. (27) recovers the
well-known textbook expression for BEC in 3D [50, 51, 62], while Eq. (26) is
in an agreement with the previously obtained results in 2D [63, 18]. In this
limit, the fraction of the polaritons in the ground state becomes macroscopic
in the 3D case as the $T$ decreases, however, in the 2D case, this fraction is
infinitesimally small for any finite temperature.

###  5.2 Entropy

The density matrix (9) allows obtaining the entropy, $S_{N}=-{\rm
Tr}(\hat{\rho}_{N}\ln\hat{\rho}_{N})$, in 2D

$S_{N}=\ln Z_{N}+G_{2D}\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}\frac{1}{n^{2}}$
(28)

and in 3D

$S_{N}=\ln
Z_{N}+\frac{3}{2}G_{3D}\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}\frac{1}{n^{5/2}}.$
(29)

The entropy per polariton decreases as the total number of polaritons grows
(Fig. 1). This is because an increase in the number of polaritons leads to an
increase in the ratio of polaritons that are in the ground state.

Figure 2:  The second-order coherence function for the polaritons in the
ground state in 2D (a) and in 3D (c), given that there are exactly $N$
polaritons distributed in the system as a whole. The spectral width of the
emission from the polaritons in the ground state in 2D (b) and in 3D (d),
given that there are exactly $N$ polaritons in the system. $T$ is the
temperature of the environment, $T_{0}$ is the temperature corresponding to
$G_{2D}=1$ in 2D and $G_{3D}=1$ in 3D (see Section 4). The black dashed line
marks the approximate condition for a macroscopic occupation of the ground
state $N=G_{2D}\ln N$ in 2D (see Eq. (26)) and $N=G_{3D}\zeta(3/2)$ in 3D (see
Eq. (27)).

###  5.3 First-order coherence function

The first-order coherence function is defined as [24, 64]

$g^{(1)}_{N}({\bf r})=\frac{\langle\hat{\psi}^{\dagger}({\bf
0})\hat{\psi}({\bf r})\rangle_{N}}{\sqrt{\langle\hat{\psi}^{\dagger}({\bf
0})\hat{\psi}({\bf 0})\rangle}_{N}\sqrt{\langle\hat{\psi}^{\dagger}({\bf
r})\hat{\psi}({\bf r})\rangle}_{N}},$ (30)

where $\hat{\psi}({\bf r})$ is a wave function with the following plane-wave
expansion:

$\hat{\psi}({\bf r})=\frac{1}{\sqrt{V}}\sum_{\bf k}\hat{a}_{\bf k}e^{i{\bf
k}{\bf r}}.$ (31)

From Eq. (9) it follows that
$\langle{\hat{a}_{\bf{k}}^{\dagger}{{\hat{a}}_{{\bf{k^{\prime}}}}}}\rangle_{N}=\left\langle{{{\hat{n}}_{\bf{k}}}}\right\rangle_{N}{\delta_{{\bf{kk^{\prime}}}}}$.
Thus, we obtain

$g^{(1)}_{N}({\bf r})=\frac{\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}}{N}+\\\
+\frac{G_{2D}}{N}\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}\frac{1}{n}e^{-r^{2}T/4\alpha
n}$ (32)

in 2D and

$g^{(1)}_{N}({\bf r})=\frac{\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}}{N}+\\\
+\frac{G_{3D}}{N}\sum_{n=1}^{N}\frac{Z_{N-n}}{Z_{N}}\frac{1}{n^{3/2}}e^{-r^{2}T/4\alpha
n}$ (33)

in 3D. The definition (30) implies that ${g^{(1)}_{N}({\bf 0})=1}$.

Above the condensation threshold, the term $\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}/N$ does not vanish, therefore, a non-equilibrium transition to
BEC is always accompanied by the buildup of long-range spatial correlations.
We note, that this result does not contradict the Mermin–Wagner theorem, that
was initially formulated for the magnetic systems [65]. In 3D, the
condensation threshold can be reached both for finite $V$ and infinite $V$,
that results in the buildup of the long-range correlations. However, in 2D,
the buildup of long-range correlations is possible only for finite $V$, but
not for infinite $V$. Indeed, according to (26) the term $\langle\hat{n}_{{\bf
k}={\bf 0}}\rangle_{N}/N$ vanishes in the limit $N,V\to\infty$,
$N/V\to\rm{const}$. Thus, for infinite sized 2D Bose gas, the condensation
threshold cannot be reached and the long-range correlations cannot arise.

Below the condensation threshold, $g^{(1)}_{N}({\bf r})=0$ as $r\to+\infty$.
This decay of second-order coherence function is exponential in the limit
$N\to\infty$, $V\to\infty$, $N/V\to{\rm const}$ with the coherence length

$l_{2D}=\sqrt{\frac{\alpha}{T}\frac{1}{\ln\left[(V/N)(T/4\pi\alpha)\right]}}$
(34)

in 2D and

$l_{3D}=\sqrt{\frac{\alpha}{T}\frac{1}{\ln\left[(V/N)(T/4\pi\alpha)^{3/2}\right]}}$
(35)

in 3D. Indeed, from Eqs. (32) and (33) by using expressions (20)–(21),
replacing the sum over $n$ by the integral and applying the Laplace’s method,
we obtain

$g^{(1)}_{N}({\bf r})\approx\frac{G_{2D}}{N}\sqrt{\frac{2\pi
l_{2D}}{r}}e^{-{r}/{l_{2D}}}$ (36)

for $r\gg l_{2D}$ in 2D and

$g^{(1)}_{N}({\bf
r})\approx\frac{G_{3D}}{N}\sqrt{4\pi\ln\left(\frac{G_{3D}}{N}\right)}\frac{l_{3D}}{r}e^{-{r}/{l_{3D}}}$
(37)

for $r\gg l_{3D}$ in 3D. The asymptotic behaviors described by Eq. (36) and
(37) are in an agreement with [62, 66].

###  5.4 Second-order coherence function

Using Eq. (12), we obtain the second-order coherence function of the ground
state

$(g^{(2)}_{{\bf k}={\bf 0}})_{N}=\frac{\langle\hat{n}_{{\bf k}={\bf
0}}\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}-\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}}{\langle\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}^{2}},$ (38)

where $\langle\hat{n}_{{\bf k}={\bf 0}}\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}=\sum_{n=1}^{N}Z_{N-n}(2n-1)/Z_{N}$. Eqs. (38) and (24)–(23)
allow us to consider the second-order coherence function in some limiting
cases. In the deep quantum regime ($T/T_{0}\ll 1$ in Fig. 2), $(g^{(2)}_{{\bf
k}={\bf 0}}(0))_{N}\approx 1-1/N$, that corresponds to a Fock state. When BEC
is not formed ($1\ll N\ll T/T_{0}$ in Fig. 2), $(g^{(2)}_{{\bf k}={\bf
0}}(0))_{N}\approx 2$, that corresponds to an incoherent state. When BEC is
formed ($1\ll T/T_{0}\ll N$ in Fig. 2), $(g^{(2)}_{{\bf k}={\bf
0}}(0))_{N}\approx 1+G_{2D}/N$ in 2D and $(g^{(2)}_{{\bf k}={\bf
0}}(0))_{N}\approx 1+4G_{3D}/N^{3/2}$ in 3D, that corresponds to a system
state tends to coherent state.

Thus, we conclude, that for a fixed total number of polaritons, the
macroscopic occupation of the ground state is accompanied by a buildup of
coherence (see Fig. 1 and Fig. 2). When $P_{N}$ (see Eq. (8)) is well
localized over some $N_{0}\gg 1$, the analysis presented above for exactly
$N_{0}$ polaritons distributed in the system as a whole is a good physical
approximation. However, in many cases, the distribution of $P_{N}$ and,
consequently, pumping scheme significantly alters the coherent properties of
the BEC. This information can be accessed by the consideration of the kinetics
of the condensation (see Section 6).

###  5.5 Spectral width

The stationary emission spectrum, $I_{\bf k}(\omega)$, of the polaritons with
the wave vector $\bf k$ is determined by the correlation function
$\langle\hat{a}_{\bf k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle$ [64, 39]

$I_{\bf k}(\omega)=\lim_{t\to+\infty}{\rm
Re}\int_{0}^{+\infty}\langle\hat{a}_{\bf k}^{\dagger}(t)\hat{a}_{\bf
k}(t+\tau)\rangle e^{-i\omega\tau}d\tau.$ (39)

According to the quantum regression theorem [64, 39], $\langle\hat{a}_{\bf
k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle={\rm Tr}\left(\hat{a}_{\bf
k}\hat{\rho}(\tau)\hat{a}_{\bf k}^{\dagger}\right)$ at $t\to\infty$ and

$\frac{d}{d\tau}{\rm Tr}\left(\hat{a}_{\bf k}\hat{\rho}(\tau)\hat{a}_{\bf
k}^{\dagger}\right)={\rm Tr}\left(\hat{a}_{\bf k}L_{\rm
pump}(\hat{\rho}(\tau)\hat{a}_{\bf k}^{\dagger})\right)\\\ +{\rm
Tr}\left(\hat{a}_{\bf k}L_{\rm diss}(\hat{\rho}(\tau)\hat{a}_{\bf
k}^{\dagger})\right)+{\rm Tr}\left(\hat{a}_{\bf k}L_{\rm
therm}(\hat{\rho}(\tau)\hat{a}_{\bf k}^{\dagger})\right),$ (40)

where at the initial moment in time $\hat{\rho}(0)$ is defined by the
stationary solution (52)–(53). To obtain $\langle\hat{a}_{\bf
k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle$, we solve the Eq. (40) in the
matrix form and then find the average $\langle\hat{a}_{\bf
k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle$. In the fast thermalization
limit, we may neglect $\hat{a}_{\bf k}L_{\rm
pump}(\hat{\rho}(\tau)\hat{a}_{\bf k}^{\dagger})$ and $\hat{a}_{\bf k}L_{\rm
diss}(\hat{\rho}(\tau)\hat{a}_{\bf k}^{\dagger})$ because these terms are
small compared to $\hat{a}_{\bf k}L_{\rm therm}(\hat{\rho}(\tau)\hat{a}_{\bf
k}^{\dagger})$. As a result, we obtain the following equation

$\hat{a}_{\bf k}\frac{d\hat{\rho}(\tau)}{d\tau}\hat{a}_{\bf
k}^{\dagger}\approx\hat{a}_{\bf k}L_{\rm therm}(\hat{\rho}(\tau))\hat{a}_{\bf
k}^{\dagger}+L_{\rm therm}^{\prime}(\hat{\rho}(\tau)),$ (41)

where

$L_{\rm therm}^{\prime}(\hat{\rho})=\hat{a}_{\bf k}\hat{a}_{\bf
k}^{\dagger}\sum_{{\bf q}\neq{\bf k}}\Gamma_{\bf qk}\hat{a}_{\bf
q}\hat{\rho}\hat{a}_{\bf q}^{\dagger}\\\ +\frac{1}{2}\sum_{{\bf q}\neq{\bf
k}}\Gamma_{\bf kq}\hat{a}_{\bf k}\hat{\rho}\hat{a}_{\bf q}\hat{a}_{\bf
q}^{\dagger}\hat{a}_{\bf k}^{\dagger}+\frac{1}{2}\sum_{{\bf q}\neq{\bf
k}}\Gamma_{\bf qk}\hat{a}_{\bf k}\hat{\rho}\hat{a}_{\bf
q}^{\dagger}\hat{a}_{\bf q}\hat{a}_{\bf k}^{\dagger}.$ (42)

To solve this equation, we assume that $L_{\rm therm}$ in the right-hand side
of Eq. (41) leads to a decay with the characteristic thermalization rate
$\Gamma$. We also assume that $L^{\prime}_{\rm therm}$ leads to a decay with a
rate much lower than $\Gamma$, but much higher than the dissipation rate and
the pumping rate. After obtaining the solution of Eq. (41), we explicitly
check these assumptions. Employing these assumptions, we search for a matrix
$\rho(\tau)$ satisfying Eq. (41) in the form (8). As a result, we obtain

$\langle\hat{a}_{\bf k}^{\dagger}(t)\hat{a}_{\bf
k}(t+\tau)\rangle_{N}=\langle\hat{n}_{\bf k}(t)\rangle_{N}e^{-i\omega_{\bf
k}\tau-(\Delta\omega_{\bf k})_{N}\tau/2}$ (43)

as $t\to+\infty$, where

$(\Delta\omega_{\bf k})_{N}=\sum_{{\bf q}\neq{\bf k}}\Gamma_{\bf kq}-\\\
-\sum_{{\bf q}\neq{\bf k}}(\Gamma_{\bf qk}-\Gamma_{\bf
kq})\frac{\langle(1+\hat{n}_{\bf k})\hat{n}_{\bf q}\rangle_{N-1}}{\langle
1+\hat{n}_{\bf k}\rangle_{N-1}}.$ (44)

We consider in more detail the linewidth of the ground state,
$(\Delta\omega_{{\bf k}={\bf 0}})_{N}$. For simplicity, we assume that the
thermalization rates do not depend on the states, i.e., $\Gamma_{\bf
q0}=\Gamma$. Note, that the thermalization rates obey the
Kubo–Martin–Schwinger relation, $\Gamma_{\bf 0q}=\Gamma\exp(-(\omega_{\bf
q}-\omega_{\bf 0})/T)$. From Eq. (12), (15) and (44) we obtain

$(\Delta\omega_{{\bf k}={\bf
0}})_{N}=\Gamma\frac{NZ_{N}-\sum_{n=0}^{N-1}Z_{n}}{\sum_{n=0}^{N-1}Z_{n}},$
(45)

When $T$ is higher than the condensation temperature, the linewidth of the
ground state emission is almost independent of the total number of polaritons
(Fig. 2). But, when $T$ is lower than the condensation temperature, an
increase in $N$ leads to the narrowing of the linewidth (Fig. 2). In the
letter case ($N\gg G_{2D}$ in 2D and $N\gg G_{3D}$ in 3D), we use Eqs. (24),
(22) and (23) to obtain

$(\Delta\omega_{{\bf k}={\bf 0}})_{N}\approx\Gamma G_{2D}\frac{\ln N}{N}$ (46)

in 2D and

$(\Delta\omega_{{\bf k}={\bf 0}})_{N}\approx\Gamma G_{3D}\frac{\zeta(3/2)}{N}$
(47)

in 3D. Thus, in 2D, the linewidth of the ground state deviates from the
Schawlow–Townes law, whereas, in 3D, the linewidth follows the Schawlow–Townes
law. The influence of the dimensionality on the Schawlow–Townes narrowing is
unexpected effect. As far as we know, this effect reported here for the first
time. It is interesting to note, that the linewidth of the emitted light for
the ground state (46) and (47) is proportional to the fraction of the
polaritons in the excited states (see Eqs. (26) and (27)).

From Eqs. (39) and (43) it follows that $\langle\hat{a}_{\bf
k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle=\sum_{N}P_{N}\langle\hat{a}_{\bf
k}^{\dagger}(t)\hat{a}_{\bf k}(t+\tau)\rangle_{N}$, therefore, the linewidth
depends on the distribution of the probabilities $P_{N}$. Eq. (44) is a good
physical approximation for a linewidth if the distribution $P_{N}$ has a sharp
peak over some $N_{0}\gg 1$. In this case, the linewidth of the emission is
Lorentzian. In other cases, the kinetics of polaritons may alter the shape of
the spectrum of the emitted light.

##  6 Kinetics of the condensation

Generally, there are two conditions for the formation of the condensate. The
first condition is the macroscopic occupation of the ground state. This
condition means that $\langle\hat{N}\rangle$ and $T$ are in the condensation
region in the $\\{N,T\\}$ plane (see Fig. 1). The second condition is the
buildup of coherence in the ground state. This condition means that the
second-order coherence function reaches $1$. Below, we show that these two
conditions differ in non-equilibrium BEC in the fast thermalization limit.

Due to dissipation and pumping processes, all the probabilities $P_{N}$ differ
from zero and may depend on time. To obtain $P_{N}$ we substitute the
thermalized density matrix (8) into Eq. (5). As a result, we obtain

$\frac{\partial P_{0}(t)}{\partial
t}=\frac{d_{0}}{Z_{1}}P_{1}(t)-\frac{p_{0}}{Z_{0}}P_{0}(t),$ (48)
$\frac{\partial P_{N}(t)}{\partial t}=\frac{d_{N}}{Z_{N+1}}P_{N+1}(t)-\\\
-\frac{d_{N-1}+p_{N}}{Z_{N}}P_{N}(t)+\frac{p_{N-1}}{Z_{N-1}}P_{N-1}(t)$ (49)

for $N>0$. The coefficients $p_{N}$ and $d_{N}$ are defined by

${d_{N}}=Z_{N+1}\sum_{\bf
k}\left({{\gamma_{\bf{k}}}+{\kappa_{\bf{k}}}}\right)\langle\hat{n}_{\bf
k}\rangle_{N+1},$ (50) ${p_{N}}=Z_{N+1}\sum_{\bf
k}{\kappa_{\bf{k}}}e^{(\omega_{\bf k}-\omega_{\bf 0})/T}\langle\hat{n}_{\bf
k}\rangle_{N+1}.$ (51)

If dissipation and pumping rates do not change in time, then the stationary
solution of Eqs. (48)–(49) can be obtained in the form of the recurrence
relation

$P_{1}=\frac{Z_{1}}{Z_{0}}\frac{p_{0}}{d_{0}}P_{0},$ (52)
$P_{N+1}=\frac{Z_{N+1}}{Z_{N}}\frac{d_{N-1}+p_{N}}{d_{N}}P_{N}-\frac{Z_{N+1}}{Z_{N-1}}\frac{p_{N-1}}{d_{N}}P_{N-1}$
(53)

for $N\geq 1$. The probability $P_{0}$ in (52) can be obtained from the
normalization condition $\sum_{N=0}^{+\infty}P_{N}=1$.

Once we know the probabilities $P_{N}$, we can find the reduced density
matrix, $\hat{\rho}_{\bf k}$, of any state with the wave vector $\bf k$. This
can be done by tracing out all the states in the Eq. (8) except for $\bf k$
(see Appendix C)

$\hat{\rho}_{\bf k}=\sum_{n=0}^{+\infty}|n\rangle\langle n|e^{-n(\omega_{\bf
k}-\omega_{\bf 0})/T}\times\\\
\times\sum_{N=n}^{+\infty}P_{N}\frac{Z_{N-n}-Z_{N-n-1}e^{-(\omega_{\bf
k}-\omega_{\bf 0})/T}}{Z_{N}}$ (54)

where state $|n\rangle$ satisfies $\hat{n}_{\bf k}|n\rangle=n|n\rangle$ and we
formally put $Z_{-1}=0$. The Eq. (54) is in agreement with the previously
obtained results [67, 55, 68]. The reduced density matrix of the ground state
has a particularly simple form

$\hat{\rho}_{{\bf k}={\bf 0}}=\sum_{n=0}^{+\infty}|n\rangle\langle
n|\sum_{N=n}^{+\infty}P_{N}\frac{Z_{N-n}-Z_{N-n-1}}{Z_{N}}$ (55)

Since Eq. (55) is valid for any set of probabilities $P_{N}$, then the
positive semi-definiteness of $\hat{\rho}_{{\bf k}={\bf 0}}$ requires
$Z_{N+1}\geq Z_{N}$ for any $N$. In Appendix A, we proof that $Z_{N}$ is
monotonic. We emphasize, that monotonicity of $Z_{N}$ is a direct consequence
of the ground state separation, that we did in the Eq. (14) (see also Eqs.
(69), (71) and (72)).

From Eq. (8) it follows, that both the particle number in the ground state
$\langle\hat{n}_{{\bf k}={\bf 0}}\rangle$ and the entropy $S$ can be expressed
through the probabilities $P_{N}$:

$\langle\hat{n}_{{\bf k}={\bf
0}}\rangle=\sum_{N=0}^{+\infty}P_{N}\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N},$ (56)
$S=\sum_{N=0}^{+\infty}P_{N}S_{N}-\sum_{N=0}^{+\infty}P_{N}\ln P_{N},$ (57)

where $\langle\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}$ is defined by Eq. (25) and
$S_{N}$ is defined by Eq. (28) in 2D and by Eq. (29) in 3D. This means that
the total average number of polaritons is the sum over the invariant subspaces
weighted with the probabilities of the system being in the corresponding
subspace.

Unlike the average number of polaritons in the ground state, $g^{(1)}({\bf
r})$ and $g^{(2)}_{{\bf k}={\bf 0}}(0)$ are not the weighted sum of
$g^{(1)}_{N}({\bf r})$ and $(g^{(2)}_{{\bf k}={\bf 0}}(0))_{N}$ over invariant
subspaces, but instead

$g^{(1)}({\bf
r})=\frac{\sum_{N=0}^{+\infty}P_{N}\langle\hat{\psi}^{\dagger}({\bf
0})\hat{\psi}({\bf
r})\rangle_{N}}{\sqrt{\sum_{N=0}^{+\infty}P_{N}\langle\hat{\psi}^{\dagger}({\bf
0})\hat{\psi}({\bf
0})\rangle_{N}}\sqrt{\sum_{N=0}^{+\infty}P_{N}\langle\hat{\psi}^{\dagger}({\bf
r})\hat{\psi}({\bf r})\rangle_{N}}},$ (58) $g^{(2)}_{{\bf k}={\bf
0}}(0)=\frac{\sum_{N=0}^{+\infty}P_{N}\langle\hat{n}_{{\bf k}={\bf
0}}\hat{n}_{{\bf k}={\bf 0}}-\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}}{\left(\sum_{N=0}^{\infty}P_{N}\langle\hat{n}_{{\bf k}={\bf
0}}\rangle_{N}\right)^{2}}.$ (59)

Therefore, the coherence of the condensate is determined not only by
$\langle\hat{N}\rangle$ and $T$, but also strongly depends on the particular
distribution of the $P_{N}$.

Below we consider in more detail the buildup of the coherence of stationary
non-equilibrium BEC at $T=0$ in the fast thermalization limit.

#### Coherence buildup at $T=0$.

As it was shown in Section 5, $\langle\hat{n}_{{\bf k}={\bf 0}}\rangle_{N}=N$,
$\langle\hat{n}_{{\bf k}\neq{\bf 0}}\rangle_{N}=0$ and $(g^{(2)}_{{\bf k}={\bf
0}})_{N}=1-N^{-1}$ at $T=0$ both in 2D and in 3D. This means that the
polaritons in ${\bf k}={\bf 0}$ in each of the invariant subspace are in the
Fock state. Nevertheless, this does not imply a buildup of the coherence
because the dissipation and the pumping processes should be taken into
account.

The coherence of the polaritons in the ground state is completely determined
by the reduced density matrix of the ground state, $\hat{\rho}_{{\bf k}={\bf
0}}$. From Eq. (55), it follows that this reduced density matrix at ${T=0}$
takes the form

$\hat{\rho}_{{\bf k}={\bf 0}}=\sum_{n=0}^{+\infty}|n\rangle\langle n|P_{n}.$
(60)

The stationary solution (52)–(53) for probabilities $P_{N}$, in this case,
reduces to

$P_{N}=\frac{1}{N!}\left(\frac{\kappa_{{\bf k}={\bf 0}}}{\gamma_{{\bf k}={\bf
0}}+\kappa_{{\bf k}={\bf 0}}}\right)^{N}\\\ \times\left(1-\frac{\kappa_{{\bf
k}={\bf 0}}}{\gamma_{{\bf k}={\bf 0}}+\kappa_{{\bf k}={\bf
0}}}\right)^{-\left(1+\kappa_{1}/\kappa_{{\bf k}={\bf 0}}\right)}\\\ \times
f_{N}\left(1+\frac{\kappa_{1}}{\kappa_{{\bf k}={\bf 0}}}\right),$ (61)

where $\gamma_{{\bf k}={\bf 0}}$ and $\kappa_{{\bf k}={\bf 0}}$ are the
dissipation and the pumping rates for ${\bf k}={\bf 0}$, and $\kappa_{1}$ is
the total pumping rate of the all excited states

$\kappa_{1}=\sum_{{\bf k}\neq{\bf 0}}\kappa_{\bf k}.$ (62)

The function $f_{N}$ is determined by

$f_{0}(x)=1,$ (63) $f_{N}(x)=x(x+1)...(x+N-1)$ (64)

for $N>0$.

The reduced density matrix (60) and probabilities (61) allow us to find the
occupation of the ground state, first- and second-order correlation functions
of the polaritons

$\langle\hat{n}_{{\bf k}={\bf 0}}\rangle=\frac{\kappa_{{\bf k}={\bf
0}}+\kappa_{1}}{\gamma_{{\bf k}={\bf 0}}},$ (65) $g^{(1)}({\bf r})=1,$ (66)
$g^{(2)}_{{\bf k}={\bf 0}}(0)=\frac{2\kappa_{{\bf k}={\bf
0}}+\kappa_{1}}{\kappa_{{\bf k}={\bf 0}}+\kappa_{1}}.$ (67)

This means that regardless of which polariton state is pumped, all the
polaritons are thermalized to the ground state and long-range spatial
correlations are established. But different pumping schemes result in
different statistics of the condensate. Depending on the pumping scheme, the
second-order coherence function of the ground state may be between $1$ and
$2$.

When only the ground state is pumped incoherently, the coherence buildup does
not occur at $T=0$. This is because the transitions from the ground state to
the excited states are suppressed and the ground state is effectively
decoupled from the excited states at $T=0$. As a result, the ground state
inherits the statistics of the incoherent pumping.

The more the excited states are pumped compared to the ground state the more
the polariton condensate is coherent. In the limiting case $\kappa_{{\bf
k}={\bf 0}}=0$, $\kappa_{1}\neq 0$, the polaritons are perfectly coherent at
$T=0$ and the reduced density matrix (60) becomes

$\hat{\rho}_{{\bf k}={\bf 0}}=\sum_{n=0}^{+\infty}|n\rangle\langle
n|\frac{1}{n!}\left(\frac{\kappa_{1}}{\gamma_{{\bf k}={\bf
0}}}\right)^{n}e^{-\kappa_{1}/\gamma_{{\bf k}={\bf 0}}}.$ (68)

This state is a randomly phased coherent state. The thresholdless condensation
into a perfect coherent state at $T=0$, in this case, is a consequence of the
fast thermalization limit. Indeed, the critical number of particles necessary
for the coherence buildup is determined by the ratio between thermalization
rate and the dissipation rate [33]. We consider the fast thermalization limit,
therefore, this critical number of particles is almost zero. Thus, we have a
thresholdless transition to a perfect coherent state (68).

##  7 Conclusion

In conclusion, we developed an analytical framework to describe non-
equilibrium ideal Bose–Einstein condensates (BEC) in the fast thermalization
limit. In this limit, we obtained an expression for the full density matrix of
non-equilibrium ideal BEC.

We showed that a macroscopic occupation of the ground state occurs under the
same conditions for the full thermal equilibrium case as for the non-
equilibrium case in the fast thermalization limit. The macroscopic occupation
of the ground state is always accompanied by the formation of long-range
first-order spatial correlations. Moreover, for a given system, only the
temperature and the average number of the particles determine whether the
ground state is macroscopically occupied. In contrast, the buildup of second-
order coherence in BEC is not fully determined by the occupation of the ground
state, but also depends on the pumping scheme. For instance, at zero
temperature, if excited states are pumped, then a coherence of the ground
state is formed, but if only the ground state is pumped, then the coherence
buildup of ground state does not occur.

At fixed temperature, an increase in the number of polaritons leads to an
increase in the total entropy, but the entropy per polariton decreases. This
is because the polaritons in the ground state have low entropy and an increase
in the number of particles leads to an increase of the fraction of the
particles in the ground state.

Above the condensation threshold, the linewidth of the ground state emission
narrows. At low temperatures and large numbers of polaritons, this narrowing
follows the Schawlow–Townes law in 3D, but, in 2D, the linewidth decreases
slower than predicted by the Schawlow–Townes law. As the temperature
decreases, the condensate line width decreases faster in 3D than in 2D.

In this work, we considered only 2D and 3D polaritons with quadratic
dispersion and non-degenerate states enumerated by the wave vectors. However,
the developed theory can be straightforwardly adopted to other dispersions,
such as dispersion, with a bottleneck, an inflection point and degenerate
states with definite wave vectors. Moreover, the description of the kinetics
is tolerant to the more realistic situation, than was considered in this work,
for instance, when radiative lifetime is bound by Hopfield expressions.

The developed theory is a theory for the ideal gas. However, this theory can
be applied to non-ideal gases, with the Hamiltonian $\hat{H}+\hat{U}$, where
$\hat{U}$ is non-ideal contribution to the energy. If the operator $\hat{U}$
commutes with the operator of the total number of polaritons $\sum_{\bf
k}\hat{n}_{\bf k}$, then the developed theory can be directly applied by
replacing the partition function for the ideal gas $Z_{N}$ with
$Z_{N}\exp(-\langle\hat{U}\rangle_{N}/T)$. However, if the non-ideal part of
the energy operator does not commute with the total number of polaritons, then
the generalization of the presented theory might be limited to the first-order
corrections.

## Acknowledgments

This work was supported by the Russian Science Foundation (Grant No.
20-72-10145). E.S.A. and V.Yu.Sh. thank the Foundation for the Advancement of
Theoretical Physics and Mathematics “Basis”. Yu.E.L. acknowledges the Basic
Research Program at the National Research University HSE.

## References

  * Plumhof et al. [2014] J.D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R.F. Mahrt.  Room-temperature Bose–Einstein condensation of cavity exciton–polaritons in a polymer.  _Nature Materials_ , 13(3):247–252, 2014.  doi: https://doi.org/10.1038/nmat3825.
  * 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.  A room-temperature organic polariton transistor.  _Nature Photonics_ , 13(6):378–383, 2019.  doi: https://doi.org/10.1038/s41566-019-0392-8.
  * Zasedatelev et al. [2021] A.V. Zasedatelev, A.V. Baranikov, D. Sannikov, D. Urbonas, F. Scafirimuto, V.Yu. Shishkov, E.S. Andrianov, Y.E. Lozovik, U. Scherf, T. Stöferle, R.F. Mahrt, and P.G. Lagoudakis.  Single-photon nonlinearity at room temperature.  _Nature_ , 597:493–497, 2021.  doi: https://doi.org/10.1038/s41586-021-03866-9.
  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen.  The road towards polaritonic devices.  _Nature Materials_ , 15(10):1061–1073, 2016.  doi: https://doi.org/10.1038/nmat4668.
  * Keeling and Kéna-Cohen [2020] J. Keeling and S. Kéna-Cohen.  Bose–Einstein condensation of exciton-polaritons in organic microcavities.  _Annual Review of Physical Chemistry_ , 71:435–459, 2020\.  doi: https://doi.org/10.1146/annurev-physchem-010920-102509.
  * Deng et al. [2003] H. Deng, G. Weihs, D. Snoke, J. Bloch, and Y. Yamamoto.  Polariton lasing vs. photon lasing in a semiconductor microcavity.  _Proceedings of the National Academy of Sciences_ , 100(26):15318–15323, 2003.  doi: https://doi.org/10.1073/pnas.2634328100.
  * Kasprzak et al. [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J.M. Keeling, F.M. Marchetti, M.H. Szymańska, R. André, J.L. Staehli, and V. Savona.  Bose–Einstein condensation of exciton polaritons.  _Nature_ , 443(7110):409–414, 2006.  doi: https://doi.org/10.1038/nature05131.
  * 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.  Exciton–polariton condensates.  _Nature Physics_ , 10(11):803–813, 2014.  doi: https://doi.org/10.1038/nphys3143.
  * Wertz et al. [2010] E. Wertz, L. Ferrier, D.D. Solnyshkov, R. Johne, D. Sanvitto, A. Lemaître, I. Sagnes, R. Grousson, A.V. Kavokin, P. Senellart, and G. Malpuech.  Spontaneous formation and optical manipulation of extended polariton condensates.  _Nature Physics_ , 6(11):860–864, 2010.  doi: https://doi.org/10.1038/nphys1750.
  * Balili et al. [2007] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West.  Bose–Einstein condensation of microcavity polaritons in a trap.  _Science_ , 316(5827):1007–1010, 2007.  doi: https://doi.org/10.1126/science.1140990.
  * 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, and A.G. Truscott.  Single-shot condensation of exciton polaritons and the hole burning effect.  _Nature Communications_ , 9(1):1–9, 2018.  doi: https://doi.org/10.1038/s41467-018-05349-4.
  * 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.  Bose–Einstein condensation of long-lifetime polaritons in thermal equilibrium.  _Physical Review Letters_ , 118(1):016602, 2017\.  doi: https://doi.org/10.1103/PhysRevLett.118.016602.
  * Deng et al. [2010] H. Deng, H. Haug, and Y. Yamamoto.  Exciton-polariton Bose-Einstein condensation.  _Reviews of Modern Physics_ , 82(2):1489, 2010\.  doi: https://doi.org/10.1103/RevModPhys.82.1489.
  * Klaas et al. [2018] M. Klaas, E. Schlottmann, H. Flayac, F.P. Laussy, F. Gericke, M. Schmidt, M. V. Helversen, J. Beyer, S. Brodbeck, H. Suchomel, et al.  Photon-number-resolved measurement of an exciton-polariton condensate.  _Physical Review Letters_ , 121(4):047401, 2018\.  doi: https://doi.org/10.1103/PhysRevLett.121.047401.
  * 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(6):4250, 1996.  doi: https://doi.org/10.1103/PhysRevA.53.4250.
  * Wei et al. [2019] M. Wei, S.K. Rajendran, H. Ohadi, L. Tropf, M.C. Gather, G.A. Turnbull, and I.D.W. Samuel.  Low-threshold polariton lasing in a highly disordered conjugated polymer.  _Optica_ , 6(9):1124–1129, 2019.  doi: https://doi.org/10.1364/OPTICA.6.001124.
  * Weill et al. [2019] R. Weill, A. Bekker, B. Levit, and B. Fischer.  Bose–Einstein condensation of photons in an erbium–ytterbium co-doped fiber cavity.  _Nature Communications_ , 10(1):1–6, 2019.  doi: https://doi.org/10.1038/s41467-019-08527-0.
  * 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ä.  Bose–Einstein condensation in a plasmonic lattice.  _Nature Physics_ , 14(7):739–744, 2018.  doi: https://doi.org/10.1038/s41567-018-0109-9.
  * 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ä.  Sub-picosecond thermalization dynamics in condensation of strongly coupled lattice plasmons.  _Nature Communications_ , 11(1):1–12, 2020.  doi: https://doi.org/10.1038/s41467-020-16906-1.
  * Malpuech et al. [2002] G. Malpuech, A. Di Carlo, A. Kavokin, J.J. Baumberg, M. Zamfirescu, and P. Lugli.  Room-temperature polariton lasers based on gan microcavities.  _Applied Physics Letters_ , 81(3):412–414, 2002\.  doi: https://doi.org/10.1063/1.1494126.
  * Banyai and Gartner [2002] L. Banyai and P. Gartner.  Real-time Bose–Einstein condensation in a finite volume with a discrete spectrum.  _Physical Review Letters_ , 88(21):210404, 2002\.  doi: https://doi.org/10.1103/PhysRevLett.88.210404.
  * Cao et al. [2004] H.T. Cao, T.D. Doan, D.B.T. Thoai, and H. Haug.  Condensation kinetics of cavity polaritons interacting with a thermal phonon bath.  _Physical Review B_ , 69(24):245325, 2004.  doi: https://doi.org/10.1103/PhysRevB.69.245325.
  * Doan et al. [2008] T.D. Doan, H. Thien Cao, D.B.T. Thoai, and H. Haug.  Coherence of condensed microcavity polaritons calculated within boltzmann-master equations.  _Physical Review B_ , 78(20):205306, 2008.  doi: https://doi.org/10.1103/PhysRevB.78.205306.
  * Tassone et al. [1997] F. Tassone, C. Piermarocchi, V. Savona, A. Quattropani, and P. Schwendimann.  Bottleneck effects in the relaxation and photoluminescence of microcavity polaritons.  _Physical Review B_ , 56(12):7554, 1997.  doi: https://doi.org/10.1103/PhysRevB.56.7554.
  * Kirton and Keeling [2013] P. Kirton and J. Keeling.  Nonequilibrium model of photon condensation.  _Physical Review Letters_ , 111(10):100404, 2013\.  doi: https://doi.org/10.1103/PhysRevLett.111.100404.
  * Kirton and Keeling [2015] P. Kirton and J. Keeling.  Thermalization and breakdown of thermalization in photon condensates.  _Physical Review A_ , 91(3):033826, 2015.  doi: https://doi.org/10.1103/PhysRevA.91.033826.
  * Strashko et al. [2018] A. Strashko, P. Kirton, and J. Keeling.  Organic polariton lasing and the weak to strong coupling crossover.  _Physical Review Letters_ , 121(19):193601, 2018\.  doi: https://doi.org/10.1103/PhysRevLett.121.193601.
  * Arnardottir et al. [2020] K.B. Arnardottir, A.J. Moilanen, A. Strashko, P. Törmä, and J. Keeling.  Multimode organic polariton lasing.  _Physical Review Letters_ , 125(23):233603, 2020\.  doi: https://doi.org/10.1103/PhysRevLett.125.233603.
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti.  Quantum fluids of light.  _Reviews of Modern Physics_ , 85(1):299, 2013\.  doi: https://doi.org/10.1103/RevModPhys.85.299.
  * Kavokin et al. [2017] A. Kavokin, J.J. Baumberg, G. Malpuech, and F.P. Laussy.  _Microcavities_.  Oxford university press, 2017.
  * Laussy [2012] F.P. Laussy.  Exciton polaritons in microcavities: New frontiers.  In Daniele Sanvitto and Vladislav Timofeev, editors, _Springer Series in solid–state sciences_ , volume 172, chapter 1, pages 1–42. Springer Science & Business Media, 2012.
  * Laussy et al. [2004a] F.P. Laussy, G. Malpuech, and A. Kavokin.  Spontaneous coherence buildup in a polariton laser.  _Physica Status Solidi (c)_ , 1(6):1339–1350, 2004a.  doi: https://doi.org/10.1002/pssc.200304064.
  * Laussy et al. [2004b] F.P. Laussy, G. Malpuech, A. Kavokin, and P. Bigenwald.  Spontaneous coherence buildup in a polariton laser.  _Physical Review Letters_ , 93:016402, Jun 2004b.  doi: https://doi.org/10.1103/PhysRevLett.93.016402.
  * Shishkov et al. [2022] V.Yu. Shishkov, E.S. Andrianov, A.V. Zasedatelev, P.G. Lagoudakis, and Y.E. Lozovik.  Exact analytical solution for the density matrix of a nonequilibrium polariton Bose–Einstein condensate.  _Physical Review Letters_ , 128:065301, Feb 2022.  doi: https://doi.org/10.1103/PhysRevLett.128.065301.
  * Lagoudakis et al. [2003] P.G. Lagoudakis, M.D. Martin, J.J. Baumberg, A. Qarry, E. Cohen, and L.N. Pfeiffer.  Electron-polariton scattering in semiconductor microcavities.  _Physical Review Letters_ , 90(20):206401, 2003\.  doi: https://doi.org/10.1103/PhysRevLett.90.206401.
  * Maragkou et al. [2010] M. Maragkou, A.J.D. Grundy, T. Ostatnickỳ, and P.G. Lagoudakis.  Longitudinal optical phonon assisted polariton laser.  _Applied Physics Letters_ , 97(11):111110, 2010\.  doi: https://doi.org/10.1063/1.3488012.
  * Coles et al. [2011] D. M. Coles, P. Michetti, C. Clark, W. C. Tsoi, A. M. Adawi, J.-S. Kim, and D. G. Lidzey.  Vibrationally assisted polariton-relaxation processes in strongly coupled organic-semiconductor microcavities.  _Advanced Functional Materials_ , 21(19):3691–3696, 2011.  doi: https://doi.org/10.1002/adfm.201100756.
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione.  _The theory of open quantum systems_.  Oxford University Press on Demand, 2002.
  * Litinskaya et al. [2004] M. Litinskaya, P. Reineker, and V.M. Agranovich.  Fast polariton relaxation in strongly coupled organic microcavities.  _Journal of Luminescence_ , 110(4):364–372, 2004\.  doi: https://doi.org/10.1016/j.jlumin.2004.08.033.
  * Mazza et al. [2009] L. Mazza, L. Fontanesi, and G.C. La Rocca.  Organic-based microcavities with vibronic progressions: Photoluminescence.  _Physical Review B_ , 80(23):235314, 2009.  doi: https://doi.org/10.1103/PhysRevB.80.235314.
  * Bittner and Silva [2012] E.R. Bittner and C. Silva.  Estimating the conditions for polariton condensation in organic thin-film microcavities.  _The Journal of Chemical Physics_ , 136(3):034510, 2012.  doi: https://doi.org/10.1063/1.3678015.
  * Ćwik et al. [2014] J.A. Ćwik, S. Reja, P.B. Littlewood, and J. Keeling.  Polariton condensation with saturable molecules dressed by vibrational modes.  _Europhysics Letters_ , 105(4):47009, 2014.  doi: https://doi.org/10.1209/0295-5075/105/47009.
  * Somaschi et al. [2011] N. Somaschi, L. Mouchliadis, D. Coles, I.E. Perakis, D.G. Lidzey, P.G. Lagoudakis, and P.G. Savvidis.  Ultrafast polariton population build-up mediated by molecular phonons in organic microcavities.  _Applied Physics Letters_ , 99(14):209, 2011.  doi: https://doi.org/10.1063/1.3645633.
  * Ramezani et al. [2018] M. Ramezani, Q. Le-Van, A. Halpin, and J.G. Rivas.  Nonlinear emission of molecular ensembles strongly coupled to plasmonic lattices with structural imperfections.  _Physical Review Letters_ , 121(24):243904, 2018\.  doi: https://doi.org/10.1103/PhysRevLett.121.243904.
  * Savvidis et al. [2000] P.G. Savvidis, J.J. Baumberg, R.M. Stevenson, M.S. Skolnick, D.M. Whittaker, and J.S. Roberts.  Angle-resonant stimulated polariton amplifier.  _Physical Review Letters_ , 84(7):1547, 2000.  doi: https://doi.org/10.1103/PhysRevLett.84.1547.
  * Kosloff [2013] R. Kosloff.  Quantum thermodynamics: A dynamical viewpoint.  _Entropy_ , 15(6):2100–2128, 2013.  doi: https://doi.org/10.3390/e15062100.
  * Shishkov et al. [2018] V.Yu. Shishkov, E.S. Andrianov, A.A. Pukhov, A.P. Vinogradov, and A.A. Lisyansky.  Zeroth law of thermodynamics for thermalized open quantum systems having constants of motion.  _Physical Review E_ , 98(2):022132, 2018.  doi: https://doi.org/10.1103/PhysRevE.98.022132.
  * Berman et al. [2008] O.L. Berman, Yu.E. Lozovik, and D.W. Snoke.  Theory of Bose–Einstein condensation and superfluidity of two-dimensional polaritons in an in-plane harmonic potential.  _Physical Review B_ , 77(15):155317, 2008.  doi: https://doi.org/10.1103/PhysRevB.77.155317.
  * Toda et al. [1983] M. Toda, R. Kubo, and N. Saito.  _Statistical Physics I: Equilibrium Statistical Mechanics_.  Springer Series in Solid-State Sciences №30. Springer Berlin Heidelberg, 1983.
  * Landau and Lifshitz [2013] L.D. Landau and E.M. Lifshitz.  _Course of theoretical physics_.  Elsevier, 2013.
  * Fraser [1951a] A.R. Fraser.  XV. The condensation of a perfect Bose–Einstein gas.–II.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 42(325):165–175, 1951a.  doi: https://doi.org/10.1080/14786445108560984.
  * Borrmann and Franke [1993] P. Borrmann and G. Franke.  Recursion formulas for quantum statistical partition functions.  _The Journal of Chemical Physics_ , 98(3):2484, 1993.  doi: https://doi.org/10.1063/1.464180.
  * Brosens et al. [1997] F. Brosens, J.T. Devreese, and L.F. Lemmens.  Thermodynamics of coupled identical oscillators within the path-integral formalism.  _Physical Review E_ , 55(1):227, 1997.  doi: https://doi.org/10.1103/PhysRevE.55.227.
  * Weiss and Wilkens [1997] C. Weiss and M. Wilkens.  Particle number counting statistics in ideal Bose gases.  _Optics Express_ , 1(10):272–283, 1997.  doi: https://doi.org/10.1364/OE.1.000272.
  * Chase et al. [1999] K.C. Chase, A.Z. Mekjian, and L. Zamick.  Canonical and microcanonical ensemble approaches to Bose–Einstein condensation: The thermodynamics of particles in harmonic traps.  _The European Physical Journal B-Condensed Matter and Complex Systems_ , 8(2):281–285, 1999.  doi: https://doi.org/10.1007/s100510050691.
  * Kocharovsky et al. [2006] V.V. Kocharovsky, V.V. Kocharovsky, M. Holthaus, C.H.R. Ooi, A. Svidzinsky, W. Ketterle, and M.O. Scully.  Fluctuations in ideal and interacting Bose–Einstein condensates: From the laser phase transition analogy to squeezed states and bogoliubov quasiparticles.  _Advances in Atomic, Molecular, and Optical Physics_ , 53:291–411, 2006.  doi: https://doi.org/10.1016/S1049-250X(06)53010-1.
  * Fraser [1951b] A.R. Fraser.  XIV. The condensation of a perfect Bose–Einstein gas.–I.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 42(325):156–164, 1951b.  doi: https://doi.org/10.1080/14786445108560983.
  * Schubert [1946] G. Schubert.  Zur Bose–Statistik.  _Zeitschrift für Naturforschung A_ , 1(3):113–120, 1946.  doi: https://doi.org/10.1515/zna-1946-0301.
  * Schubert [1947] G. Schubert.  Zur Bose-Statistik (Nachtrag).  _Zeitschrift für Naturforschung A_ , 2(5):250–251, 1947.  doi: https://doi.org/10.1515/zna-1947-0503.
  * Holthaus and Kalinowski [1999] M. Holthaus and E. Kalinowski.  The saddle-point method for condensed Bose gases.  _Annals of Physics_ , 276(2):321–360, 1999.  doi: https://doi.org/10.1006/aphy.1999.5950.
  * Pitaevskii and Stringari [2003] L. Pitaevskii and S. Stringari.  _Bose–Einstein Condensation_.  International series of monographs on physics 116 Oxford science publications. Clarendon Press, 2003.
  * Ketterle and Van Druten [1996] W. Ketterle and N.J. Van Druten.  Bose–Einstein condensation of a finite number of particles trapped in one or three dimensions.  _Physical Review A_ , 54(1):656, 1996.  doi: https://doi.org/10.1103/PhysRevA.54.656.
  * Scully and Zubairy [1999] M.O. Scully and M.S. Zubairy.  Quantum optics, 1999.
  * Mermin and Wagner [1966] N.D. Mermin and H. Wagner.  Absence of ferromagnetism or antiferromagnetism in one-or two-dimensional isotropic heisenberg models.  _Physical Review Letters_ , 17(22):1133, 1966\.  doi: https://doi.org/10.1103/PhysRevLett.17.1133.
  * Mouchliadis and Ivanov [2008] L. Mouchliadis and A.L. Ivanov.  First-order spatial coherence of excitons in planar nanostructures: A k-filtering effect.  _Physical Review B_ , 78(3):033306, 2008.  doi: https://doi.org/10.1103/PhysRevB.78.033306.
  * Wilkens and Weiss [1997] M. Wilkens and C. Weiss.  Particle number fluctuations in an ideal Bose gas.  _Journal of Modern Optics_ , 44(10):1801–1814, 1997.  doi: https://doi.org/10.1080/09500349708231847.
  * Mossel and Caux [2012] J. Mossel and J.-S. Caux.  Exact time evolution of space-and time-dependent correlation functions after an interaction quench in the one-dimensional Bose gas.  _New Journal of Physics_ , 14(7):075006, 2012\.  doi: https://doi.org/10.1088/1367-2630/14/7/075006.
  * Comtet [2012] L. Comtet.  _Advanced Combinatorics: The art of finite and infinite expansions_.  Springer Science & Business Media, 2012.

##  Appendix A Some properties of the partition function $Z_{N}$

In this Appendix, we present some general properties of the partition function
$Z_{N}$. To derive these properties, we introduce the generating function of
$Z_{N}$: $\sum_{N=0}^{+\infty}x^{N}Z_{N}$, $0<x<1$. Eq. (10) allows us to
represent the generating function in different forms

$\sum_{N=0}^{+\infty}x^{N}Z_{N}=\sum_{N=0}^{+\infty}x^{N}\sum_{{\rm
config}:\sum_{\bf k}n_{\bf k}=N}e^{-\sum_{\bf k}n_{\bf k}(\omega_{\bf
k}-\omega_{\bf 0})/{T}}=\sum_{N=0}^{+\infty}\sum_{{\rm config}:\sum_{\bf
k}n_{\bf k}=N}e^{-\sum_{\bf k}n_{\bf k}((\omega_{\bf k}-\omega_{\bf
0})/{T}-\ln x)}=\\\ =\prod_{\bf k}\sum_{n_{\bf k}=0}^{+\infty}e^{-n_{\bf
k}((\omega_{\bf k}-\omega_{\bf 0})/{T}-\ln x)}=\prod_{\bf
k}\frac{1}{1-xe^{-(\omega_{\bf k}-\omega_{\bf 0})/{T}}}=\exp\left\\{-\sum_{\bf
k}\ln\left[1-xe^{-(\omega_{\bf k}-\omega_{\bf 0})/T}\right]\right\\}.$ (69)

The last two forms are especially useful for the considerations below.

###  A.1 Proof of Eq. (12)

From Eq. (69) it follows that

$\sum_{N=0}^{+\infty}x^{N}\frac{\partial Z_{N}}{\partial\omega_{\bf
q}}=-\frac{1}{T}\frac{xe^{-(\omega_{\bf q}-\omega_{\bf
0})/T}}{1-xe^{-(\omega_{\bf q}-\omega_{\bf 0})/T}}\prod_{\bf
k}\frac{1}{1-xe^{-(\omega_{\bf k}-\omega_{\bf 0})/{T}}}=\\\
=-\frac{1}{T}\frac{xe^{-(\omega_{\bf q}-\omega_{\bf
0})/T}}{1-xe^{-(\omega_{\bf q}-\omega_{\bf
0})/T}}\sum_{N=0}^{+\infty}x^{N}Z_{N}=\sum_{N=0}^{+\infty}x^{N}\frac{1}{T}\left(Z_{N}-\sum_{n=0}^{N}Z_{N-n}e^{n(\omega_{\bf
q}-\omega_{\bf 0})/T}\right).$ (70)

Eq. (12) directly follows from this equation.

###  A.2 Monotonicity of $Z_{N}$

To prove the monotonicity of $Z_{N}$ we use Eq. (69) and consider the
expression

$\sum_{N=0}^{+\infty}x^{N}Z_{N}=\frac{1}{1-x}\sum^{+\infty}_{N=0}x^{N}F_{N},$
(71)

where we introduced

$\sum^{+\infty}_{N=0}x^{N}F_{N}=\prod_{{\bf k}\neq{\bf
0}}\frac{1}{1-xe^{-(\omega_{\bf k}-\omega_{\bf 0})/{T}}}.$ (72)

Although finding $F_{N}$ is hard in the general case, the positivity of
$F_{N}$ follows directly from Eq. (72). Moreover, from Eq. (71) it follows
that $Z_{N}=\sum^{N}_{n=0}F_{n}$. Thus, $Z_{N}$ rises monotonously with $N$.

###  A.3 The limit of $Z_{N}$ as $N\to+\infty$

Below, we assume that the limit of $Z_{N}$ as $N\to+\infty$ exists. In this
case, we introduce the notation $Z_{\infty}=\lim_{N\to+\infty}Z_{N}$

$Z_{\infty}=\lim_{N\to+\infty}\left[\sum_{n=1}^{N}(Z_{n}-Z_{n-1})+Z_{0}\right]=\lim_{x\to
1-0}\lim_{N\to+\infty}\left[\sum_{n=1}^{N}x^{n}(Z_{n}-Z_{n-1})+xZ_{0}\right]=\\\
=\lim_{x\to
1-0}\lim_{N\to+\infty}\left[\sum_{n=0}^{N}x^{n}Z_{n}-x\sum_{n=0}^{N-1}x^{n}Z_{n}\right]=\lim_{x\to
1-0}(1-x)\sum_{N=0}^{+\infty}x^{N}Z_{N}.$ (73)

Thus, we express $Z_{\infty}$ through the generating function of $Z_{N}$.

###  A.4 An exact expression for $Z_{N}$ in 2D and in 3D for quadratic
dispersion

To find an exact expression for $Z_{N}$ we use its generating function (see
Eq. (69))

$\sum_{N=0}^{+\infty}x^{N}Z_{N}=\exp\left\\{-\sum_{\bf
k}\ln\left[1-xe^{-(\omega_{\bf k}-\omega_{\bf 0})/T}\right]\right\\}$ (74)

Below we consider the quadratic dispersion, $\omega_{\bf k}=\omega_{\bf
0}+\alpha{\bf k}^{2}$, in cases of 2D and 3D.

#### Polaritons in 2D.

In 2D, the density of states is $\nu(\omega)=V/4\pi\alpha$ and Eq. (74) takes
the form

$\sum_{N=0}^{+\infty}x^{N}Z_{N}=\frac{\exp\left[G_{2D}\leavevmode\nobreak\
{\rm Ln}_{2}(x)\right]}{1-x},$ (75)

where $G_{2D}=VT/4\pi\alpha$ is the number of states in the energy range $T$
above ${\bf k}={\bf 0}$ and ${\rm Ln}_{2}(x)$ is the polylogarithm of order
$2$.

In Eq. (75), we use the expansion of the exponent in the series of Bell
polynomials [69] and the expansion of $(1-x)^{-1}$ in a Taylor series, then we
collect the coefficients with equal powers of $x$ on the left and right sides
and obtain

$Z_{N}=\sum_{n=0}^{N}\frac{1}{n!}B_{n}\left(\frac{1!\leavevmode\nobreak\
G_{2D}}{1^{2}},\frac{2!\leavevmode\nobreak\
G_{2D}}{2^{2}},\leavevmode\nobreak\ ...\leavevmode\nobreak\
,\frac{n!\leavevmode\nobreak\ G_{2D}}{n^{2}}\right),$ (76)

where $B_{n}$ is the $n$-th complete exponential Bell polynomial [69] and
$B_{0}=1$. Thus, in 2D, $Z_{N}$ is fully determined by a dimensionless
parameter $G_{2D}$ and rises monotonically from $Z_{0}=1$ to
$Z_{\infty}=\exp(G_{2D}\zeta(2))$, where $\zeta(2)$ is the value of the zeta
function at $2$.

#### Polaritons in 3D.

In 3D, the density of states is
$\nu(\omega)=(V/4\pi^{2}\alpha)\sqrt{(\omega-\omega_{\bf 0})/\alpha}$ and Eq.
(74) takes the form

$\sum_{N=0}^{+\infty}x^{N}Z_{N}=\frac{\exp\left[G_{3D}\leavevmode\nobreak\
{\rm Ln}_{5/2}(x)\right]}{1-x},$ (77)

where $G_{3D}=V(T/4\pi\alpha)^{3/2}$ is the number of states in the energy
range $T$ above ${\bf k}={\bf 0}$ and ${\rm Ln}_{5/2}(x)$ is the polylogarithm
of the order $5/2$.

From Eq. (77) it follows that

$Z_{N}=\sum_{n=0}^{N}\frac{1}{n!}B_{n}\left(\frac{1!\leavevmode\nobreak\
G_{3D}}{1^{5/2}},\frac{2!\leavevmode\nobreak\
G_{3D}}{2^{5/2}},\leavevmode\nobreak\ ...\leavevmode\nobreak\
,\frac{n!\leavevmode\nobreak\ G_{3D}}{n^{5/2}}\right),$ (78)

Thus, in 3D, $Z_{N}$ is fully determined by the dimensionless parameter
$G_{3D}$ and rises monotonically from $Z_{0}=1$ to
$Z_{\infty}=\exp(G_{3D}\zeta(5/2))$, where $\zeta(5/2)$ is the value of the
zeta function at $5/2$.

#### The calculation $Z_{N}$.

We note that the direct application of Eq. (76) and (78) is not the fastest
way to calculate $Z_{N}$. This is due to the great computational difficulty of
the complete Bell functions for $N\gg 1$. The optimal strategy to obtain
$Z_{N}$ is to use the recurrence relation, discussed in the main text.

##  Appendix B Proof of Eqs. (26) and (27)

According to Eq. (25), when $N\gg 1$,

$\frac{N-\langle\hat{n}_{{\bf k}={\bf
0}}\rangle}{N}=\frac{\sum_{n=0}^{N-1}(Z_{N}-Z_{n})}{NZ_{N}}\approx\frac{\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})}{NZ_{\infty}}.$
(79)

We consider $\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})$ in more detail.

In 2D, we have

$\lim_{N\to+\infty}\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})=\lim_{N\to+\infty}\lim_{x\to
1-0}\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})x^{n}=\\\ =\lim_{x\to
1-0}\frac{\exp(G_{2D}{\rm Ln}_{2}(1))-\exp(G_{2D}{\rm
Ln}_{2}(x))}{1-x}=Z_{\infty}G_{2D}\lim_{x\to 1-0}{\rm
Ln}_{1}(x)=Z_{\infty}G_{2D}\lim_{N\to+\infty}\ln(N).$ (80)

Therefore, when $N\gg 1$

$\frac{N-\langle\hat{n}_{{\bf k}={\bf 0}}\rangle}{N}\approx\frac{G_{2D}}{N}\ln
N$ (81)

Eq. (26) directly follows from Eq. (81).

In 3D, we have

$\lim_{N\to+\infty}\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})=\lim_{N\to+\infty}\lim_{x\to
1-0}\sum_{n=0}^{N-1}(Z_{\infty}-Z_{n})x^{n}=\\\ =\lim_{x\to
1-0}\frac{\exp(G_{3D}{\rm Ln}_{5/2}(1))-\exp(G_{3D}{\rm
Ln}_{5/2}(x))}{1-x}=Z_{\infty}G_{3D}\lim_{x\to 1-0}{\rm
Ln}_{3/2}(x)=Z_{\infty}G_{3D}\zeta(3/2).$ (82)

Therefore, when $N\gg 1$,

$\frac{N-\langle\hat{n}_{{\bf k}={\bf
0}}\rangle}{N}\approx\frac{G_{3D}}{N}\zeta(3/2).$ (83)

Eq. (27) directly follows from Eq. (83).

##  Appendix C Derivation of the Eq. (54)

To derive Eq. (54), we denote the trace over all the states, except for the
state $\bf k$, as ${\rm Tr}^{\prime}_{\bf k}$. This notation allow us to write
the reduced density matrix $\hat{\rho}_{\bf k}$ of the state with the wave
vector $\bf k$ in terms of the full density matrix $\hat{\rho}$ in the simple
form

$\hat{\rho}_{\bf k}={\rm Tr}^{\prime}_{\bf k}\hat{\rho}$ (84)

Using the Eq. (8), we obtain

$\hat{\rho}_{\bf k}=\sum_{N=0}^{+\infty}P_{N}{\rm Tr}^{\prime}_{\bf
k}\hat{\rho}_{N}=\sum_{N=0}^{+\infty}P_{N}\frac{1}{Z_{N}}\sum_{n=0}^{N}|n\rangle\langle
n|\sum_{{\rm config}:\sum_{{\bf q}\neq{\bf k}}n_{\bf q}=N-n}e^{-n(\omega_{\bf
k}-\omega_{\bf 0})/{T}-\sum_{{\bf q}\neq{\bf k}}n_{\bf q}(\omega_{\bf
q}-\omega_{\bf 0})/{T}}=\\\
=\sum_{N=0}^{+\infty}P_{N}\frac{1}{Z_{N}}\sum_{n=0}^{N}|n\rangle\langle
n|e^{-n(\omega_{\bf k}-\omega_{\bf 0})/{T}}z_{N-n}({\bf
k})=\sum_{n=0}^{+\infty}|n\rangle\langle n|e^{-n(\omega_{\bf k}-\omega_{\bf
0})/{T}}\sum_{N=n}^{+\infty}P_{N}\frac{z_{N-n}({\bf k})}{Z_{N}},$ (85)

where $|n\rangle$ is the state, that satisfies $\hat{n}_{\bf
k}|n\rangle=n|n\rangle$, and $z_{N}({\bf k})$ is the partition function of all
the states but the state with wave vector $\bf k$

$z_{N}({\bf k})=\sum_{{\rm config}:\sum_{{\bf q}\neq{\bf k}}n_{\bf
q}=N}e^{-\sum_{{\bf q}\neq{\bf k}}n_{\bf q}(\omega_{\bf q}-\omega_{\bf
0})/{T}}$ (86)

By analogy with Eq. (69), we obtain the generating function for $z_{N}({\bf
k})$

$\sum_{N=0}^{+\infty}x^{N}z_{N}({\bf k})=\prod_{{\bf q}\neq{\bf
k}}\frac{1}{1-xe^{-(\omega_{\bf q}-\omega_{\bf 0})/{T}}}$ (87)

Comparing Eqs. (69) and (87), we conclude $\sum_{N=0}^{+\infty}x^{N}z_{N}({\bf
k})=(1-xe^{-(\omega_{\bf k}-\omega_{\bf
0})/{T}})\sum_{N=0}^{+\infty}x^{N}Z_{N}$. Therefore,

$z_{N}({\bf k})=Z_{N}-e^{-(\omega_{\bf k}-\omega_{\bf 0})/{T}}Z_{N-1}$ (88)

where we formally put $Z_{-1}=0$. Substitution of Eq. (88) into Eq. (85) leads
to Eq. (54).