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

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.

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 [1,6,7].
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 phononand 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 pop-ulation 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 ρ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 timedependent 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 nonequilibrium 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.

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 polari-ton systems.Scattering by electrons [36], highenergy 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 [1,7,13] We consider a continuous system of polariton states with the frequency ω k for a wave vector k in the general case of D dimensional space.We suppose that the polaritons with wave vector k of the lower polariton branch are described by bosonic creation â † k and annihilation âk operators [31,32] which obey the commutation relation âk , â † q = δ(k − q).In this case, the Hamiltonian of the polaritons takes the form where the state with k = 0 is the ground state.
Hereafter, we assume that h = 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 ρ 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] (2) where γ k is the dissipation rate of the state with wave vector k.The effective incoherent pumping of the lower polaritons can be described by the Lindblad superoperator where κ k is the pumping rate of the state with wave vector k.The Lindblad operator L pump leads to the following dynamics of the average polariton population in the state with the wave vector k [35]: (d âk /dt) pump = 0. Therefore, L 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 polaritonpolariton scattering [14,46].However, regardless of the mechanism, the thermalization can be described through the Lindblad superoperator [31,32] where Γ kq is the transition rate from the state with the wave vector k to the state with the wave vector q.The thermalization rates Γ kq obey the Kubo-Martin-Schwinger relation Γ kq /Γ qk = exp ((ω k − ω q ) /T ), 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 ρ has the form 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, Γ 0k (1+ â † 0 â0 ) κ k , γ k .In such a case, at times t γ −1 k , κ −1 k , only thermalization affects the system dynamics.As a result, the density matrix obeys the approximate differential equation After the time , the system reaches its quasistationary state determined by The relaxation operator L therm (ρ), defined by the expression (4), conserves the total number of lower polaritons [31].Thus, in the limit ) → 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 therm (ρ) conserves the total number of lower polaritons [31], the operator of total number of lower polaritons k â † k â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 k n 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: 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.ρN is the thermalized density matrix of the states with exactly N polaritons, which forms an invariant subspace [35].This means that ρN can be represented as where Z N is the partition function of the states with exactly N polaritons, and Rconfig is a diagonal density matrix corresponding to the configuration of the polaritons {n k }.This density matrix obeys â † k âk Rconfig = n k Rconfig , k n 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 nondiagonal elements while finding spectral linewidth in Section 5.
The partition function Z N is defined by the normalization condition Tr( Rconfig ) = 1, thus When there are no polaritons in the system, Eq. ( 9) implies ρN=0 = Z −1 0 |0 0|.Therefore, The density matrix (8) allows to calculate all the single-time averages.In particular, it allows calculating the averages nm The derivative of Z N with respect to ω k is (see Appendix A) (13) Note that ∂Z N /∂ω 0 must be calculated as lim k→0 ∂Z N /∂ω k .

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 k by the integral over frequencies ω: The density of states ν(ω) is 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 k = 0 is non-degenerate.The separation of the state k = 0 is a standard step in the description of continuous systems [31,50,51].
Below we consider the quadratic dispersion, ω k = ω 0 + αk 2 .This dispersion is characteristic for polaritons of the lower branch in the vicinity of k = 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, where nk N stands for Tr(n k ρN ).According to Eqs. ( 12) and (13), After integration over k we obtain the recurrence relations with G 2D = V T /4πα in 2D and with G 3D = V (T /4πα) 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 k = 0 in the corresponding dimensions.
The exact expressions for Z N are obtained in Appendix A trough Darwin-Fowler method [52,[58][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 T /T 0 in the figures).From the Eqs.(18) and (19) we obtain in 2D and 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 T /T 0 N in the figures).From the Eqs.(80) and (82) we obtain with with 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 1 in 2D and G 3D 1 (T /T 0 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 5 Condensation of exactly N polaritons 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.

Population of the states
From Eq. (12) it follows that the average number of polaritons in the state with the wave vector k can be calculated as 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 → ∞, from Eq. (25) we obtain the fraction of polaritons in the ground state (see Appendix B) in 2D and in 3D.
In the limit N, V → ∞, N/V → 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 [18,63].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.

Entropy
The density matrix (9) allows obtaining the entropy, S N = −Tr(ρ N ln ρN ), in 2D and in 3D 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.

First-order coherence function
The first-order coherence function is defined as [24,64] g (1) where ψ(r) is a wave function with the following plane-wave expansion: From Eq. (9) it follows that â † k âk N = nk N δ kk .Thus, we obtain in 2D and g (1) (33) in 3D.The definition (30) implies that g (1) Above the condensation threshold, the term nk=0 N /N does not vanish, therefore, a nonequilibrium 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 nk=0 N /N vanishes in the limit N, V → ∞, N/V → const.Thus, for infinite sized 2D Bose gas, the condensation threshold cannot be reached and the longrange correlations cannot arise.
Below the condensation threshold, g N (r) = 0 as r → +∞.This decay of second-order coherence function is exponential in the limit N → ∞, V → ∞, N/V → const with the coherence length in 2D and 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) for r l 2D in 2D and for r l 3D in 3D.The asymptotic behaviors described by Eq. ( 36) and (37) are in an agreement with [62,66].

Second-order coherence function
Using Eq. (12), we obtain the second-order coherence function of the ground state (g where nk=0 nk=0 N = N n=1 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 1 in Fig. 2), (g k=0 (0)) N ≈ 1 − 1/N , that corresponds to a Fock state.When BEC is not formed (1 N T /T 0 in Fig. 2), (g k=0 (0)) N ≈ 2, that corresponds to an incoherent state.When BEC is formed (1 T /T 0 N in Fig. 2), (g 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 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).

Spectral width
The stationary emission spectrum, I k (ω), of the polaritons with the wave vector k is determined by the correlation function â † k (t)â k (t + τ ) [39, 64] According to the quantum regression theorem [39,64] where at the initial moment in time ρ(0) is defined by the stationary solution (52)- (53).To obtain â † k (t)â k (t + τ ) , we solve the Eq.(40) in the matrix form and then find the average â † k (t)â k (t + τ ) .In the fast thermalization limit, we may neglect âk L pump (ρ(τ )â † k ) and âk L diss (ρ(τ )â † k ) because these terms are small compared to âk L therm (ρ(τ )â † k ).As a result, we obtain the following equation where To solve this equation, we assume that L therm in the right-hand side of Eq. ( 41) leads to a decay with the characteristic thermalization rate Γ.We also assume that L therm leads to a decay with a rate much lower than Γ, 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 ρ(τ ) satisfying Eq. (41) in the form (8). As a result, we obtain (43) as t → +∞, where We consider in more detail the linewidth of the ground state, (∆ω k=0 ) N .For simplicity, we assume that the thermalization rates do not depend on the states, i.e., Γ q0 = Γ.Note, that the thermalization rates obey the Kubo-Martin-Schwinger relation, Γ 0q = Γ exp(−(ω q − ω 0 )/T ).From Eq. ( 12), ( 15) and (44) we obtain 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 G 2D in 2D and N G 3D in 3D), we use Eqs.(24), ( 22) and (23) to obtain in 2D and 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) , 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 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.

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 N 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 for N > 0. The coefficients p N and d N are defined by 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 for N ≥ 1.The probability P 0 in (52) can be obtained from the normalization condition Once we know the probabilities P N , we can find the reduced density matrix, ρk , of any state with the wave vector k.This can be done by tracing out all the states in the Eq. ( 8) except for k (see Appendix C) where state |n satisfies nk |n = n|n and we formally put Z −1 = 0.The Eq. ( 54) is in agreement with the previously obtained results [55,67,68].
The reduced density matrix of the ground state has a particularly simple form 55) is valid for any set of probabilities P N , then the positive semi-definiteness of ρk=0 requires Z N +1 ≥ 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 nk=0 and the entropy S can be expressed through the probabilities P N : where nk=0 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) (r) and g ( k=0 (0) are not the weighted sum of g (1) N (r) and (g (2) k=0 (0)) N over invariant subspaces, but instead .
(59) Therefore, the coherence of the condensate is determined not only by N 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, nk=0 N = N , nk =0 N = 0 and (g k=0 ) N = 1 − N −1 at T = 0 both in 2D and in 3D.This means that the polaritons in k = 0 in each of the invariant subspace are in the Fock state.Nevertheless, this does not imply a buildup of the coherence 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, ρk=0 .From Eq. (55), it follows that this reduced density matrix at T = 0 takes the form The stationary solution (52)-( 53) for probabilities P N , in this case, reduces to where γ k=0 and κ k=0 are the dissipation and the pumping rates for k = 0, and κ 1 is the total pumping rate of the all excited states The function f N is determined by 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 g (1) (r) = 1, (66) 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 κ k=0 = 0, κ 1 = 0, the polaritons are perfectly coherent at T = 0 and the reduced density matrix (60) becomes 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).

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 nondegenerate 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 nonideal gases, with the Hamiltonian Ĥ + Û , where Û is non-ideal contribution to the energy.If the operator Û commutes with the operator of the total number of polaritons k nk , then the developed theory can be directly applied by replacing the partition function for the ideal gas Z N with Z N exp(− Û N /T ).However, if the nonideal 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.
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 : +∞ N =0 x N Z N , 0 < x < 1. Eq. (10) allows us to represent the generating function in different forms The last two forms are especially useful for the considerations below.

A.2 Monotonicity of Z N
To prove the monotonicity of Z N we use Eq.(69) and consider the expression where we introduced 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 = N n=0 F n .Thus, Z N rises monotonously with N .
A. 3 The limit of Z N as N → +∞ Below, we assume that the limit of Z N as N → +∞ exists.In this case, we introduce the notation x N Z N .(73) Thus, we express Z ∞ 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)) +∞ N =0 x N Z N = exp − k ln 1 − xe −(ω k −ω 0 )/T (74) Below we consider the quadratic dispersion, ω k = ω 0 + αk 2 , in cases of 2D and 3D.
Polaritons in 2D.In 2D, the density of states is ν(ω) = V /4πα and Eq.(74) takes the form where G 2D = V T /4πα is the number of states in the energy range T above k = 0 and 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 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 , where ζ( 2) is the value of the zeta function at 2.

Polaritons in 3D.
In 3D, the density of states is ν(ω) = (V /4π 2 α) (ω − ω 0 )/α and Eq.(74) takes the form x N Z N = exp G 3D Ln 5/2 (x) where G 3D = V (T /4πα) 3/2 is the number of states in the energy range T above k = 0 and Ln 5/2 (x) is the polylogarithm of the order 5/2.From Eq. (77) it follows that 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 1.The optimal strategy to obtain Z N is to use the recurrence relation, discussed in the main text.

Figure 1 :
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 ζ(3/2) in 3D (see Eq. (27)).

Figure 2 :
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 ζ(3/2) in 3D (see Eq.(27)).