Mean field study of 2D quasiparticle condensate formation in presence of strong decay

Bose-condensation in a system of 2D quasiparticles is considered in the scope of a microscopic model. Mean-field dynamical equations are derived with the help of the Schwinger-Keldysh formalism and a simple model is proposed which allows to describe key features of condensate formation in systems with various quasiparticle decay rates. By analysing stationary solutions of this equation, we obtain the phase diagram of quasiparticle gas, finding a bistability region in the parameter space of the system. Finally, as an application of our theory, we study the phase diagram of a 2D exciton-polariton system in CdTe microcavity.

Though during condensate formation the cold atom gas is out of equilibrium, the resulting condensate state is an equilibrium one, which is not the case for solid state quasiparticle systems.Despite the similar nature of the low-temperature state in these systems and in cold atomic gases, the quasiparticle condensate is different in several aspects.Firstly, due to finite lifetime, these systems need to be pumped externally, hence the condensate is in a quasi-equilibrium state which is determined by an interplay between pumping and decay processes.Moreover, these condensates are often considered in low dimensional systems and together with small masses of quasiparticles, it may change totally the relevant scales (energies, times) of the condensate and the process of its formation.One of the most attractive features is the potential ability to observe high temperature condensation.
The basic Gross-Pitaevskii equation, describ-N.A. Asriyan: naasriyan@vniia.ruing condensate in equilibrium systems has been modified in numerous ways (leading to dissipative Gross-Pitaevskii-type models) in order to describe phenomenologically non-equilibrium physics of exciton/photon/exciton-polariton condensates.To step beyond the description of the kinetic stage of condensate formation, which had been well studied [12][13][14] and to incorporate the coherent properties of condensate in the evolution equation, condensate is commonly considered to be an open quantum system subject to reservoirs, namely the pump reservoir, decay reservoir, lattice phonon reservoir, etc.Several phenomenological models were proposed in this scope.Namely, the most common is the one introduced by M. Wouters and I. Carusotto [15], which describes polariton condensate as a system, coupled to classical excitonic reservoir of density n R : (1) Here P describes incoherent pumping, R(n R ) is an amplification rate of reservoir-condensate scattering, g and g stand for intracondensate and condensatereservoir particle interaction.Parameters γ R and γ are decay rates of quasiparticles from reservoir and condensate.
Gain saturation, which is essential for describing condensate density equilibrating, may be incorporated directly into the condensate evolution equation: These appeared to be fruitful approaches, which allowed to describe spontaneous vortex lattice formation [16], pattern formation [17], as well as relaxation oscillations were considered [18].
Models of the type (1) and (2) are well-suited for describing long-lifetime systems with moderate decay rates γ ≪ γ R .Though, with the help of (1) the opposite case (γ ≫ γ R ) was also considered [19,20], where the modulational instability of the homogeneous condensate was demonstrated to be a consequence of slow reservoir relaxation.
However, the above-mentioned models share several drawbacks.They do not allow considering the normal phase-condensate transition in systems with wide range of decay rates of quasiparticles: from the "ultracold atom gas" limit with negligible decay to the "strong dissipation" case in presence of high losses from condensate compensated by high gain rates from the reservoir.For the latter case, one should consider the impact of the broadening of the condensate spectral function on the interaction with the reservoir, which introduces memory effects.
As we will show further, the key to overcoming this issue is considering the frequency-dependent gain.One way to do this is using a phenomenological non-Hermitian term first introduced by L.P. Pitaevskii [21] for superfluid He II and later adopted for describing frequency dependent gain due to the polaritonpolariton interaction [22]: ( Here, close to the threshold frequency Ω tr the gain efficiency decreases to zero.Using such a model for polariton condensate appeared to be essential for obtaining its excitation spectrum and describing its superfluid properties, as shown in [22]. Along with studying the model equations, microscopic theories were developed for exciton/exciton polariton systems [23,24] in order to derive evolution equations ab initio.This activity resulted in deriving microscopically motivated expressions for the terms of equations (1) and (2).However, a similar microscopical treatment, which could reproduce the frequencydependent gain and could be applicable in presence of strong decay, is lacking.
In this paper, we propose a microscopically motivated equation for condensate formation dynamics in 2D finite-lifetime quasiparticle system, which allows considering consistently the cases of various decay rates due to a properly described frequencydependent gain effect.Stated otherwise, the approach we use aims to be applicable to both the case of long-living particles (i.e.weakly coupled to the decay bath) and rapidly-decaying quasiparticles whose lifetime may be much smaller than reservoir evolution time scales.In the latter case the reservoir acts as a system with long memory effect providing a delayed feedback on condensate evolution, i.e. it may be considered as a non-Markovian evolution regime.Using this dynamical equation, we obtain the meanfield phase diagram for a quasiparticle system and demonstrate that finite lifetime of condensate particles may lead to formation of condensate and normal phase overlap regions on the phase diagram.Moreover, we show that above-mentioned memory effects significantly affect the condensate dynamics leading to evolution patterns different from exponential approach to equilibrium and density relaxation oscillations.As a demonstration of real-life application of the dynamical equation we derive, we will consider the phase diagram for CdTe microcavity polariton gas.
Being motivated by the discussion of spontaneous symmetry breaking in cold atom condensate by H. T. C. Stoof [25,26], we use the same theoretical framework: the Schwinger-Keldysh technique in path integral formulation.It provides direct access to the condensate order parameter, that's why it has been widely applied to quasiparticle condensates.For instance, photonic condensate in a dye-filled optical microcavity [27], exciton polariton condensates in quantum wells [23,24], parametrically pumped polariton systems [28] were also considered.We use here a single-level condensate model (e.g.assume uniform condensate in a finite-sized system) in the spirit of pioneering articles [12,14] dealing with exciton condensation.Of course, this prevents us from studying the spatial structure of the condensed state, namely from reproducing the results of [20].However, this simplification allows keeping the further discussion analytically tractable and focusing on the impact of intense condensate decay on the evolution equations.
We start from defining in Sec. 2 the model Hamiltonian for 2D quasiparticles with pair interaction.Here our goal is to develop a general theory, not considering any particular system.Introducing the necessary parameters, we derive the mean-field dynamical equations for condensate evolution in Sec. 3.This derivation is followed by a description of a possible simplified approximate expressions for equation terms.It leads to a dynamical model, which is one of the main results of the current paper.
The next Sec.4 is devoted to discussion of a meanfield phase diagram for a generic quasiparticle system.Possible phases are described followed by a stability analysis.A detailed description of condensation dynamics is presented in Sec. 5.
After a discussion on the advantages and drawbacks of the model in Sec.6, in the last Sec.7 of the article we present an application of the model to real-life system such as exciton polariton gas in CdTe/CdMgTe microcavity.

The model system
We deal with a 2D single-level condensate model of finite-lifetime quasiparticles embedded in a long-living particle reservoir.The system is treated as a 2D bose gas with contact interparticle interaction and leakage from the condensate.The model Hamiltonian for the system described is as follows ( ψq is the quasiparticle annihilation operator with wavevector q, therefore ψ0 corresponds to the condensate mode): Here g 0 = V0 2L 2 stands for contact interparticle interaction with V 0 being the interaction potential and L 2 denoting the quantization area.ε denotes the possible energy detuning of the condensate level with respect to the reservoir dispersion curve ε q→0 = 0.A decay rate Γ is introduced to describe condensate particle finite lifetime.Hereafter ℏ = 1.
The action for this system defined on the Schwinger-Keldysh contour in path integral formulation is as follows: Our goal is to integrate out the reservoir degrees of freedom to derive an effective action for the condensate.We are going to treat the reservoir in the simplest possible way as a continuously pumped quasiequilibrium system with stationary surface density n and effective temperature T = (k B β) −1 .Moreover, we assume the collision broadening for the reservoir to be negligible compared to the effective temperature.When considering the greater/lesser components of the Green's function on the Schwinger-Keldysh contour, this allows to use the following approximation (the spectral function is assumed to be a sharply peaked Lorentzian, see Fig. 1): Here f (x) = (e x − 1) −1 , ζ q = ε q +∆ε q −µ with ∆ε q taking into account the blueshift due to reservoir interparticle interaction.And µ stands for the ideal gas chemical potential of the reservoir.From the expressions above, we may combine the casual Green's function: with γ res → 0. Note that the assumptions introduced in this section (the form of the spectral function with uniform broadening of all the states as well as γ res → 0) are reasonable for the case of weak interparticle interaction in the reservoir.For real systems, this corresponds to weakly interacting gases.For instance, it is the case for exciton-polariton systems at reasonable temperatures.A rough estimate is V0n 2D kT ≪ 1.For typical excitonic densities of order n 2D ∼ 10 10 cm −2 this leads to T ≫ 2µK.
The necessary assumption of "narrow" spectral function.The red and blue lines schematically depict the frequency dependencies of the bosonic distribution function and the spectral function of a reservoir quasiparticle, respectively.

Effective action
We may integrate out reservoir degrees of freedom in order to obtain an effective action for the condensate mode only (with ψ being the corresponding field), which has the following structure: Here we introduced a self-energy term Σ to describe interaction with reservoir, extracting the time-local contribution Σ δ (τ, τ ′ ) = Σ δ δ(τ, τ ′ ) explicitly.We use the lowest order diagrammatic expressions for these terms, as presented in the Fig. 2.
In order to deal with real-time dynamics, we pass from the field ψ defined on the Keldysh contour (denoted by c in (8)) to ψ ± being the fields on its backward and forward branches.This is achieved by the standard Keldysh rotation This procedure makes the real-time action to acquire the following form: with local (Σ δ 0 ), retarded (Σ + ), advanced (Σ − 0 ) and Keldysh (Σ K 0 ) components of the self-energy term introduced.q , given by (7).
To derive the dynamical equation from this type of action, one needs to treat noise terms coupled to the ϕ field dynamics, which leads to a Langevin type equation.In the current paper, we focus on the mean-field dynamics by seeking the stationary phase "classical" solution in a form: with ϕ c (t) being a solution of the following equation: Analysing this equation is the main objective of the current work.
Neglecting the noise term imposes several limitations.Namely, we are not able to describe correctly the initial stages of condensate evolution when fluctuations dominate the dynamics.The approach presented below leads to relevant results in the vicinity of stationary points, where noise terms are less significant.Incorporating them terms into the theory is left for future investigations.

Self-energy term
As we see in the Fig. 2, the self-energy term besides an evident blueshift contribution (with N res and n res being the total reservoir occupation and its surface density correspondingly) has a term due to interparticle interaction (here τ = t − t ′ ): Note that the term Σ δ q is independent of q, therefore reservoir states are shifted by the same amount and the energy offset ε is not affected by Σ δ .In fact, by this we use the Hartree-Fock approximation for the reservoir (for contact interaction direct and exchange terms give the same contribution, that's why we have a factor of 2 on the first line in the Fig. 2).
For the retarded component, one obtains the following expression by setting q = 0 (partly following [24] by two of us): Considering γ res → 0 as well as taking advantage of the Bose-Einstein distribution property 1 + f (ω) = e ω f (ω), we obtain the following expression for the imaginary part of the retarded component: Above A dimensionless function I(ω) may be isolated as follows: with λ dB = ℏ √ 2mT being the thermal de-Broglie wavelength.Hereafter β = 1, which means all the energies are measured in units of kT .
One may evaluate I(ω) numerically, considering the quadratic dispersion relation.Moreover, asymptotic behaviour may be studied analytically (see Appendix A to find calculations for arbitrary momentum q).All the information is summarized in the Fig. 3.The dash-dotted line qualitatively demonstrates how the curve will be modified for finite-range interaction.Here for illustration µ = −0.2 is chosen.
When discussing the ω → ∞ asymptotics, we note that I(ω) = const does not vanish.This is due to the contribution of the following process: a virtual particle with q = 0 and ε = w scatters on a reservoir particle with q ≈ 0, ϵ q ≈ 0. As a result, they are both in the reservoir with momenta q 1 = −q 2 and energies The amplitude of this process does not decrease with growing ω.This is a consequence of the contact interaction model.If we consider some finite interaction radius, i.e., introduce a transfer momentum cut-off, I(ω) will tend to zero with growing ω because of growing transferred momentum q ∼ √ ω during the scattering process.This is schematically demonstrated by a dash-dotted red line in the Fig.

The exact behaviour depends, of course, on the interparticle interaction potential.
Using this numerical result for further calculations is quite involved.We will further use an approximation which captures the double-peaked shape of the curve as well as reproduces the x-intercept correctly, which appears to be important to describe condensate effective chemical potential equilibrating.Here, Λ + and γ should be treated as fitting parameters.
Of course, the model curve cannot reproduce the real one exactly.As it will be clear from Sec. 6, when searching the best fit parameters, it's worth better fitting the left peak at the cost of not reproducing the high-frequency behaviour as demonstrated in the Fig. 3.
The retarded self-energy term is analytic in the upper half plane in frequency domain due to causality, therefore its real and imaginary part are related by the Kramers-Kronig relations, and we may seek an approximation for the function using its imaginary part only.Therefore: In time domain, this function is as follows (note that time is measured in units of 1/kT ): 4 Mean-field phase diagram

Reduction to an ODE system
Using (21), we may write down the dimensionless dynamical equation as follows: The advantage provided by using an exponential kernel is the possibility to simplify the dynamical equation (22).To do that, one may consider it along with its time derivative in order to exclude the memory term.Namely, the Madelung transformation which leads to the ODE system after differentiating and excluding integral terms: To get rid of negative power terms in ρ, we may now use ρ(t) = e r(t) substitution in order to obtain the following autonomous ODE system: By introducing V = ṙ and ν = θ as new variables, we may formulate the dynamical equation as a first order ODE system: To derive initial conditions, one may set t = 0 in (23): This leads to ν(0) = ϵ + ge r and V (0) = −Γ.Note that these expressions are due to assuming the condensate and reservoir being in contact from t = 0 exactly.In real system the initial evolution stages may be more complicated, which, however, does not change the asymptotic behaviour of the system.

Stationary points
Here, we seek stationary points of the ODE system (27) by setting We readily solve these equations to obtain: For stability analysis of these stationary points, one may consider the linear expansion and obtain the corresponding eigenvalues (see Appendix B).The "upper" solution with "+" sign appears to be stable and the "lower" one -unstable.

Decaying solution
Dealing with a constrained quantity ρ ≥ 0, we should consider one more equilibrating scenario with ρ → 0. We may seek the decaying solution in a form r = 0, ṙ = κ, θ = Ω = const.With small ρ being neglected, this leads to the following characteristic equations: As derived in Appendix B, there are two eigenmodes with κ < 0 (i.e., the decaying solution is stable) whenever ρ + ρ − ≥ 0 or D < 0. Stated otherwise, the decay solution is stable if there are either no stationary points or both of them are present simultaneously.

Phase diagram
From the discussion above, we infer that there are two possible equilibrating scenarios.The one is reaching a stationary solution ρ + and the other is a decaying solution.In physical terms, the first one corresponds to condensate formation (with non-zero ϕ being the corresponding order parameter) and the second one describes the normal phase.
Using the results of the stability analysis, we may summarize them on a phase diagram presented in the Fig. 4 on (Λ + , Γ) plane (Λ + physically corresponds to condensate-reservoir interaction "strength", Γ is the decay intensity of condensate particles).The remaining parameters γ and µ are fixed.
On this figure, the condensate exists whenever ρ + > 0 (which implies D > 0).The more strict condition of these two (ρ + > 0 and D > 0) defines the condensate stability boundary.
The decaying solution is stable in the three cases listed below: 1.For regions with D < 0.Here no stationary points exists, decay is the only asymptotic scenario.This is the "Normal phase" region of the diagram, below the straight line.
2. When D > 0 but ρ − < 0 as well as ρ + < 0. There are also no physically relevant stationary points here, this is another part of the "Normal phase" region, which is the in the left bottom corner of the diagram, above the dashed line.Note that for low enough, ε this region disappears (e.g. in the Fig. 6).
3. For D > 0, ρ + > 0 and ρ − > 0. Here, both the decaying solution and one of the stationary points are stable.This is the "Bistability" region on the phase diagram bounded by D = 0 and ρ − = 0 lines (note that ρ − < ρ + , that's why it is ρ − < 0 which makes the stability criterion ρ − ρ + > 0 invalid).This is kind of an overlap of "Normal phase" and "Condensate" regions of the phase diagram.
The equations of the phase boundaries are presented near the corresponding lines.
There is kind of a triple point on the phase diagram where all the three solutions coexist with ρ ± → 0. Its position is given by and regardless of the detuning, ε it is located on the straight line by blue arrows on the dashed line in the Fig. 4).Note that for Γ → 0, only the "upper" stationary point exists.It is given by: The effective chemical potential ν of the condensate becomes equal to the one of the particle reservoir, as

Condensate Bistability
Normal phase one could expect for atomic gas of long-living particles.The condensation threshold may be identified at ρ → 0 (see the y-intercept in the Fig. 4): Since Λ + , µ and γ themselves are not independent quantities, but they depend on density and temperature, this equation may be treated as a one defining the critical effective temperature T eff c .When increasing condensate decay rate, this temperature gets shifted.A demonstration of critical temperature evaluation will be presented in Sec.7 where we map this phase diagram on the density/temperature plane.

Evolution in different regimes
In the two of the three phases described, condensate formation is possible.In the "Condensate" phase there is a single stationary point present which attracts all the ODE solutions regardless of the initial conditions as demonstrated in the Fig. 5 (a).
In contrast, in the "Normal phase" region all the solutions are attracted towards ρ = 0 with θ approaching ν ∞ , which corresponds to the slower decaying eigenmode as illustrated in In the bistability regime, an unstable stationary point appears which repels the occupation to either ρ = 0 or the stable point as illustrated in the Fig. 5  (c)-(d).Note that in the Fig. 5 (c) the upper stationary point also exists at a higher occupation.One may infer, studying the evolution in the bistable regime, that higher initial occupations are attracted to the stationary point (condensate formation takes place as in the Fig. 5 (d)) and lower ones decay to zero as presented in the Fig. 5 (c).However, it is not always the case, since for some parameters even for large ρ the line of initial condition ν = ε + gρ does not intersect the attraction basin of the stable point.

Condensate formation. Relaxation oscillations
On fig 5 (a) one may see oscillations when approaching equilibrium.Though, such type of relaxation oscillations are not a general feature of the system.
As it is shown in details in Appendix B, the eigenmode expansion of ϕ(t) close to stationary state consists of terms e κ0t , e (κ1±iΩ)t .We expect significant asymptotic density oscillations similar to the ones described in the scope of a different model in [18], when κ 0 < κ 1 (note that both are negative) in order for the oscillatory terms to dominate at late times.This regime is illustrated in the Fig. 7 (a)-(c).In contrast, for κ 0 > κ 1 , at late times condensate occupation monotonously approaches stationary value as shown in the Fig. 7 (d)-(f).
Relaxation oscillations are present for low enough ε, the typical phase diagram for this regime is presented in the Fig. 6.
Of course, the oscillations may be observable outside the hatched region of the phase diagram also.However, at late times, they are replaced by monotonous exponential approach to the stationary state.Note that oscillations are damped for higher values of Γ, which may be considered as a consequence of increasing the impact of memory terms by means of decreasing the lifetime of the particles in the condensate.

Bifurcation scheme
In order to draw several physical predictions in the bistable region, we study here the bifurcation diagram of the ODE system.It is presented schematically in the Fig. 8 for fixed Γ, µ, g and ε.It may be treated as a cross-section of the phase diagram Fig. 4 by a vertical line passing to the right of the "triple point".
We readily observe the bistable region where there are two stable branches.Along with the repelling behaviour of the lower stationary point, illustrated in the Fig. 5 (c)-(d), we expect two physical effects in this region.
Foremost, when condensate is formed in the bistability region, hysteresis is possible when changing the Λ + by varying the quasiparticle density.This is illustrated by a cycle of violet arrows in the Fig. 8.
The second prediction is related to condensate formation dynamics.One should note that initial conditions for ϕ c cannot be defined precisely, at least due to uncertainty relation.We may set up the initial distribution only.This makes condensate formation a probabilistic process.In any particular realization, with all the other parameters being the same, the system may end up either in normal or in condensed phase.
Detailed study of both these effects is beyond the scope of the mean-field analysis, fluctuations should be systematically treated.(f) Eigenmodes for this case.These plots correspond to ▲ in the Fig. 6.

Discussion
6.1 Alternative treatment.

The origin of phases
The integro-differential equation (13) itself may provide some useful qualitative understanding even without converting it to an ODE system.Namely, considering a stationary solution in a form ϕ c = √ ρe −iνt (assuming non-zero ρ and ν = const), we may derive the following pair of equations: The first of them provides a relation for the selfconsistent effective chemical potential ν of the condensate.The second one describes particle flux saturation.Obviously, we should impose a constraint gρ ≥ 0, therefore, the necessary condition for condensation is existence of a non-empty set of zeros ν * of the imaginary part ℑ[Σ + 0 (ν)], which satisfy the condition ℜ [Σ + (ν * )] ≤ ν * − ε 0 − Σ δ .This condition may be considered as another version of the condensate formation criterion described by H.T.C. Stoof in [25] when describing spontaneous symmetry breaking in cold atomic gas.
These general statements may be illustrated with the use of the model expression for the self-energy term.The equation (22) The stationary points (ν ± , ρ ± ) may be obtained from the graphical representation of the equation system (38) below in the Fig. 9.One may identify two intersection points of ℑ[Σ + (ν)] with Γ 2 in the Fig. 9 as the ones, corresponding to the two stationary points (30).By solving the first equation of the system with ν ± with respect to ρ, we find stationary occupation numbers ρ ± , which are physically relevant if ρ ≥ 0. Graphically, ρ ± > 0 when the corresponding points A ± with coordinates (ν ± , ℜ[Σ + (ν ± )]) are in the shaded region below the line ν − ε.For instance, in the Fig. 9 only one stationary point A + is in the shaded region, which corresponds to the "Condensate" phase.
This qualitative discussion provides an illustration for the claim from Sec. 3 about focusing on the left peak when fitting the ℑ[Σ + (ω)] curve, since it defines the stationary points.This is by no means a fully justified statement, since the real part of the self energy may be significantly modified even in the vicinity of ω = 0 along with changing the high-frequency behaviour of the imaginary part.It means that our model is incapable of describing fine effects due to the exact form of the interparticle interaction (which defines mostly how exactly does the ℑ[Σ(ω)] decay at high frequencies).It is suitable for qualitative description only.
However, for an arbitrary Σ + (ω), one can apply the same graphical procedure and seek solutions in the shaded region where ν − ε ≥ ℜ[Σ(ν)].This approach allows not only equilibrating the occupation ρ = |ϕ| 2 of the condensate, but also describing the spontaneous symmetry breaking and stationary phase dynamics.

Reservoir particle interaction
As one may infer from equation system (37), the condensate formation is crucially dependent on the form of ℜ[Σ + 0 (ω)].Moreover, since reservoir levels may be shifted also, the relative offset ℜ[Σ + 0 (ω)−Σ + q (ω)] with respect to reservoir states with non-zero q is relevant.Taking this into account will change the quantitative predictions of the theory.However, as demonstrated in Appendix A, the asymptotic behaviour of Σ + (q) is independent of q, the double-peaked shape remains the same as well as Σ + q (µ) = 0 for arbitrary q.That's why we hope that a thorough treatment (the one similar to what is done for 3D cold atom gas in Ref. [30]) will not to affect qualitative predictions of the model.

Late time evolution
We may now show how the model equation ( 22) is related to the driven-dissipative Gross-Pitaevski model (2) (for uniform system since we deal with a singlelevel condensate).To do that, one should consider the particle flux dependence on the condensate occupation number.As discussed above, using an ansatz ϕ = √ ρe −iνt , one may derive the equation pair (38) from (22).By expanding near the stable stationary point (ν + , ρ + ), we express the particle flux as follows: with the coefficients given by We have here the coefficients of (2) expressed in terms of the ones of (38) in the vicinity of the stable point.
For long-living particles with Γ = 0 (which leads to ν + = µ) the expression for Γ eff is given by However, note that in contrast to (2), the model approach developed here with an exponential memory kernel describes how the reservoir imposes not only the occupation but the condensate effective chemical potential also.This is due to frequency dependent gain, which is described by the frequency dependence of the memory kernel.This is crucial for describing condensate formation and its phase dynamics.That's what allows us to identify the phase boundaries.

Numerically fitting the memory kernel
In order to adopt the presented model for describing real-life systems, one needs to perform numerical integration over the polariton momenta in (17) (see Appendix A for details) and then use fitting to evaluate Λ + and γ.Remarkably, this can be done just once since the dimensionless function I + (ω), which was introduced in Sec. 3 is only dependent on the normalized density ñ = nλ 2 db .Performing the fitting for various ñ (see details at the end of Appendix A) results in approximate expressions of the form: which are reasonably accurate for ñ ∈ [0.01; 1].

Demonstration
Above, we studied the system, given by (4).Here we will demonstrate how to adopt the results to a particular quasiparticle system such as an exciton-polariton gas.
Generally, for low enough temperatures, one may consider lower polaritons with dispersion Here Ω is the Rabi splitting, m ph and m ex stand for photon and exciton masses respectively, ∆ is the photon dispersion detuning with respect to the excitonic one.It is given as follows (E g is the semiconductor gap, E b < 0 is the exciton binding energy): with D being here the microcavity width.Condensate is mainly localized at the minimum, and the reservoir particles occupy the "flat" part of the spectrum.
Therefore, we may argue that the condensate offset is given by which is negative regardless of the sign of ∆.
|ε| Condensate Reservoir E(k) When describing polariton-polariton interaction, we need to take into account Hopfield coefficients.Namely, the excitonic coefficient is given by (k is the polariton momentum) We further use some simplifications.Namely, given that excitonic mass is usually negligible compared to the photonic one (the ratio is of order 10 −4 − 10 −3 ), we note that the momentum-dependent detuning ∆ significantly exceeds Rabi splitting Ω in case for typical reservoir momenta.Given the reservoir temperature, this assumption is valid for kT ≫ mexc m ph Ω.The typical Rabi splitting for quantum well excitons is of order Ω ∼ 10 meV, this assumption is violated only for extremely low temperatures T ≪ 1 K.
That's why we adopt this assumption (∆ ≫ Ω for reservoir) to consider the reservoir as dominantly excitonic with Since interpolariton interaction is mainly due to the excitonic component, we may use a model Hamilto-nian of the following form: This leads to an additional factor of |X 0 | 2 for Σ + .Moreover, blueshifts for the condensate and the reservoir are different now, and we need to renormalize the offset energy as follows: When considering condensate polariton decay rate, one may express it as a function of cavity photon lifetime τ : with 1 − |X 0 | 2 representing the photonic component of the condensate polariton.As a real-life example, we may treat in this fashion a polariton gas in CdTe/CdMgTe microcavites [8].The necessary parameters are Ω = 26 meV, m exc = 0.69 meV, V 0 = 1.8 µeV × µm 2 .
With the help of the approach described above, we solve numerically the equations D = 0, ρ ± = 0 (which define the phase boundaries as illustrated in the Fig. 4).
In the Fig. 11 the phase diagram is presented in terms of reservoir density and temperature for several experimentally accessible detuning values and condensate lifetimes.It shows how the lines in the Fig. 4 are mapped on density-temperature plane and allows localizing all the three phases.These plots demonstrate that the bistable region is indeed more significant for high decay rates of condensate particles (lower lifetimes), its existence is due to intense condensate decay.

Conclusion
In this paper, we considered a model for describing quasiparticle condensate as an open system embedded in a quasi-equilibrium reservoir.The corresponding open-dissipative Gross-Pitaevski type equation for the condensate has an integral memory term due to the influence of the reservoir.We proposed a simplification, which allows treating the complex integrodifferential dynamical equation as an autonomous ODE system.Dealing with stationary solutions of this On the plots "N", "B" and "C" stand for "Normal phase", "Bistability" and "Condensate" respectively.Phase boundaries are dashed in the regions where the fitting procedure and the analytical approximation (44) no more reliable (ñ is outside the region [0.01; 1]).
system, we described a phase diagram predicting the existence of a bistable phase.Several dynamical effects were described, including relaxation oscillations, hysteresis, etc.
To demonstrate the real-life applicability of the model, we considered a polariton gas in CdTe microcavity deriving a phase diagram for this model and localizing the regions of the condensed/normal phases of the system.Though not claiming full agreement with experiment, we expect the form of the phase diagram, in particular the existence of the bistable region, to be a general feature of condensates in quasiparticle systems with finite lifetime.
We see the main advantage of the proposed dynamical equation being its ability to naturally describe frequency-dependent gain and describe correct equilibrium behaviour (in sense of aligning condensate and reservoir chemical potentials) in finite-lifetime limit.
We expect the approach presented here with the structure of the memory kernel proposed to be useful for describing not only the mean field stationary states, but also dynamical, statistical properties far from equilibrating.This requires incorporating fluctuations in the model.Doing so and considering coherence build-up in the condensate is a subject of future work.

A The self-energy term A.1 Analytical calculations
We start the discussion here using an expression from the main text (note that this expression is dimensionless, β = 1 is considered): Considering quadratic spectrum ε q = q 2 2m for reservoir particles and passing to continuum limit with dimensionless momentum q → qλ dB , we obtain: It's now helpful to use q + = q1+q2 2 and q − = q 1 − q 2 as integration variables: After integrating out the delta-function: Here n − = q− q− is a unitary 2D vector.We now use q δ = q + −q and the corresponding unitary vector n δ to obtain: The theta-function specifies the lower integration limit as follows: This expression may be used for evaluating I q (ω).However, asymptotic behaviour may be studied analytically.
Considering ω → ∞: The opposite limit ω → ∞ (note that q ∼ 1 is the thermal momentum magnitude): Note that these asymptotic expressions are independent of q.
A.2 Fitting the I 0 (ω) curve We have numerically integrated the expression for I 0 (ω) from (59) and performed fitting of the left peak (as described in the main text) with a model expression.The results, as well as the approximate least-squares estimates which justify analytical approximations for the Λ + (ñ) and γ(ñ) dependences are presented below in the Table 1 and Fig. 12:

B Stability analysis B.1 Stationary points. Stability
For stability analysis, we need to express equations in the form of an ODE system: We may linearise it by substitution ν = ν eq + δν, V = δV, r = r eq + δr.Leaving first order terms only: Seeking the eigenmodes of this system in a form Using expressions (30) for the stationary points, one may simplify (+ sign in the last term corresponds to the upper stationary point): with f 2 ± being expressed as: First note that whenever the lower stationary point exists (D > 0 and ρ − ≥ 0) the free term on the left-hand side of (70) is negative.It is enough to conclude that the lower solution has at least one positive real eigenvalue χ > 0, leading to instability.For the upper stationary point, all the terms are positive whenever it exists.Therefore, the cubic equation (70) has one real negative eigenvalue and two complex-conjugate eigenvalues.We may explicitly consider the real and imaginary components of χ = κ + iΩ, which leads to the following equation system: The non-zero frequency satisfying the second equation is Ω 2 = 3κ 2 + 2κ (2γ + Γ) + f 2 ± .By substitution, we get: The free term here may be expressed as follows: We may now note that all the coefficients of the cubic equation (73) are positive for + sign chosen, which leads to κ < 0 and allows to finally conclude that the upper stationary point is stable.

B.2 Decaying solutions
We start from the characteristic equations for the decaying solution, presented in the main text: After excluding ν from these equations and using a substitution κ → −γ − Γ 2 + y: Since the last term is given by subtraction of a non-negative perfect square, in order for all the roots to be less than γ + Γ 2 (which is the same as requesting all the κs to be negative), we need to impose a condition f Γ 2 + γ ≥ 0. By direct substitution, one may verify that this stability condition is equivalent to the inequality or ρ + ρ − ≥ 0. From the expression above, we see that D < 0 (when there are no stationary points at all) also satisfies the stability condition.The eigenvalues themselves are of the following form: The corresponding frequencies are:

Figure 2 :
Figure 2: Time-local term and the contribution to the selfenergy term due to interaction with reservoir.Here the filled vertex corresponds to the factor (−g0), solid lines are reservoir Green's functions G res q , given by (7).

Figure 3 :
Figure 3: Frequency dependence of the retarded self-energy term component I(ω) (red solid line) with asymptotics given.The dash-dotted line qualitatively demonstrates how the curve will be modified for finite-range interaction.Here for illustration µ = −0.2 is chosen.

Figure 4 :
Figure 4: Phase diagram for the system under consideration.The blue dot denotes the "triple" point of coexistence of all the three phases.The arrows indicate how does this point move with varying ε.Here for demonstration γ = 0.2, g = 1, µ = −0.2,ε = 0.7.For each of the black dots (a)-(d) there is an evolution graph presented in the Fig. 5.
Fig 5 (b) (see Appendix B for details).

Figure 5 :
Figure 5: Evolution illustrated for (a) "Condensate" phase; (b) "Normal phase"; (c)-(d) "Bistability" phase for two initial occupations.The letters correspond to points in the Fig. 4. The dashed line restricts the initial conditions by ν = ε + gρ as given by (28).The hollow circle indicates the initial occupation chosen.

Figure 8 :
Figure 8: Bifurcation diagram for the bistable region of the phase diagram.The stable solution branches are drawn as solid curves, for the unstable one a thick dashed curve is used.The bifurcation point is pictured in white.Arrows indicate its movement along the dashed line with Γ being changed.The double arrow on the right indicates how does the line move with ε being varied.

Figure 10 :
Figure 10: Polariton dispersion.For the lower polariton branch, the energy detuning ε of the condensate with respect to the reservoir is indicated.

Figure 11 :
Figure 11: Phase diagrams for CdTe microcavity polaritons for various photon lifetimes τ and detunings ∆.On the plots "N", "B" and "C" stand for "Normal phase", "Bistability" and "Condensate" respectively.Phase boundaries are dashed in the regions where the fitting procedure and the analytical approximation (44) no more reliable (ñ is outside the region [0.01; 1]).