Interplay of coherent and dissipative dynamics in condensates of light

Based on the Lindblad master equation approach we obtain a detailed microscopic model of photons in a dye-filled cavity, which features condensation of light. To this end we generalise a recent non-equilibrium approach of Kirton and Keeling such that the dye-mediated contribution to the photon-photon interaction in the light condensate is accessible due to an interplay of coherent and dissipative dynamics. We describe the steady-state properties of the system by analysing the resulting equations of motion of both photonic and matter degrees of freedom. In particular, we discuss the existence of two limiting cases for steady states: photon Bose-Einstein condensate and laser-like. In the former case, we determine the corresponding dimensionless photon-photon interaction strength by relying on realistic experimental data and find a good agreement with previous theoretical estimates. Furthermore, we investigate how the dimensionless interaction strength depends on the respective system parameters.


Introduction
Within the last decades open dissipative many-body quantum systems have emerged as a promising research direction for both basic research and applications. In particular, this is due to the development of exquisite technologies to coherently manipulate and control the internal and external degrees of freedom of atomic and photonic matter, as well as their interaction. A prominent example at the immediate interface of quantum optics and condensed matter physics is provided by the laser as a coherent light source which has contributed not only to our modern understanding of non-equilibrium phase transitions in general but even to many useful applications in our everyday life [1][2][3][4].
Another more modern prominent object of research is the Bose-Einstein condensate (BEC) of light, which has so far been realised in a dye-filled micro-cavity at room temperature in Bonn [5], in London [6], and quite recently also in Utrecht [7]. One of the key ingredients is the possibility of photons to acquire an effective mass by trapping them in a cavity in two dimensions -without this, the photons would just disappear according to the Planck law upon lowering the temperature. In the experiment, this is achieved via a curved-mirror cavity, which changes the dispersion relation of the photons from linear to quadratic. Along the resonator axis the frequency of the mode is quantised according to the resonance condition. Simultaneously, the curved mirrors create a harmonic trapping potential for the photons in the transverse direction. The next crucial element is given by dye molecules in the resonator, which are pumped incoherently. The multiple absorption and emission events between the cavity photons and the dye molecules lead to a thermalisation of the light [8], so the resulting photon BEC emerges from an equilibrium phase transition [9][10][11]. Photon thermalisation was also shown to be possible in much simpler but periodically driven systems, such as double quantum dots [12], or a collection of harmonic modes [13], both coupled to some environment. Recently, an elaborate theoretical proposal for a BEC of light in nanofabricated semiconductor micro-cavities has been put forward [14].
The usual atomic BEC as a thermal equilibrium phase transition occurs for temperatures below some critical value [15,16] when the resulting ground state condensate acquires macroscopic occupation, while the populations of the higher energy levels obey the Bose-Einstein distribution even before the transition. Phase transition into a macroscopically occupied mode emerges as well in a controllable way in the laser case, but − in contrast to the BEC − this transition depends on the rate of the loss and the pumping channels, which makes it a non-equilibrium paradigm. The two transitions, which a priori seem to be incompatible with each other due to their different nature, can thus be regarded as two sides of the same coin within a single non-equilibrium setup [17]. Several studies concerning the similarities and differences of condensate and lasing states and their appearance in different systems exist [18][19][20][21][22][23]. The investigation of systems and conditions under which a complex equilibrium state can be realized within a non-equilibrium setup has developed into an attractive topic, both in experiment and theory.
Apart from photonic systems, condensation effects of bosonic quasi-particles have also been observed in solid-state physics for magnons [24][25][26][27][28] and exciton polaritons [29][30][31][32][33][34]. The latter quasi-particles can be created in semiconductor micro-cavities using strong coupling between photons and particle-hole excitations [29]. A non-equilibrium BEC of polaritons has been observed in various experiments in polymers [35,36]. Surprisingly, the transition there is not always restricted to a mode with the lowest momentum [37]. The first microscopic model of a photon condensate was developed by Kirton and Keeling [38,39], which has recently been further extended by the same authors [40,41]. They considered a dye-filled cavity with multiple optical modes together with additional incoherent pump and loss channels and derived a Markovian quantum master equation of the Lindblad type [42,43]. Using an adiabatic elimination of the degrees of freedom of the dye molecules, Kirton and Keeling obtained a mean-field equation for the occupation of the cavity modes. The resulting steady state turned out to have different physical properties depending on the values of the respective system parameters. Provided that the relaxation time towards equilibrium is much shorter than the life time of the photons in the cavity, the steady state is given by a Bose-Einstein distribution, otherwise a laser-like state occurs having macroscopic occupation of a higher energy mode. Inspired by such a behaviour, a minimal two-mode laser model with a Dicke-like interaction was investigated [44]. Different phases with up to four possible and up to two stable fixed points were found, some of which have an analogy to the laser-to-condensatelike transition. However, this analogy is only quite limited due to the absence of a temperature scale in the model. Quite recently, by considering the full spatial dynamics of light [40] a rich non-equilibrium phase diagram featuring Bose-Einstein condensation, multimode condensation and lasing has been demonstrated [45]. On the other side, by using the Schwinger-Keldysh formalism a Langevin field equation describing the dynamics of photons in a dye-filled cavity was obtained [46] and later utilised to study phase fluctuations [47] and phase diffusion [48] in such systems. Moreover, a quantum Langevin model for non-equilibrium condensation of photons in planar microcavity devices was developed in [49] and recently extended to address pseudo-thermalisation in driven-dissipative non-Markovian open quantum systems [50]. A theoretical description of a photon condensate based on three-dimensional Maxwell equations, which are mapped via a paraxial approximation to a two-dimensional Schrödinger equation, was suggested as well [51]. We also note that a unified theory for excited-state, fragmented and equilibrium-like Bose condensation in pumped photonic many-body systems has recently been introduced in [52].
The theoretical modelling of dissipative condensates usually strives for a reduced description in terms of a mean-field approximation in the form of a complex-valued Gross-Pitaevskii equation, which explicitly takes into account gains and losses. It describes the system around this phase transition even in non-equilibrium [53][54][55]. Within equilibrium, a real-valued Gross-Pitaevskii equation is a standard tool to describe condensation effects [15,16,[56][57][58][59]. At the present stage, the Gross-Pitaevskiilike equation for a photon condensate can only be obtained by including a non-linear self-interaction into the model on a phenomenological level [5,49,51]. A more detailed investigation shows that this non-linear self-interaction of photons is mediated via the change of the refractive index of the dye molecules due to the mutual presence of the optical Kerr and the thermo-optical effect [60]. Due to dimensional reasons the effective photon-photon interaction strength g in two spatial dimensions corresponds to a dimensionless numberg = gm/ 2 [61], which turns out to be of the order of 10 −9 −10 −8 for the Kerr and 10 −4 for the thermo-optic effect, respectively [60]. Based on the observed momentum-and position-resolved spectra and images of the photoluminescence from thermalised and condensed dye-microcavity photons, the upper boundg 10 −3 was obtained [62]. In addition, a theoretical investigation of the influence of photon-photon interaction on the number fluctuations in a BEC of light [63] successfully explained the measurements [64] and estimated the rangeg ∼ 10 −8 − 10 −7 . Surprisingly, even a much higher value for the interactiong ∼ 10 −2 was recently measured [65].
In this paper we generalise the microscopic model of the photon BEC by Kirton and Keeling [38,39] such that the dye-mediated contribution to the photon-photon interaction strength becomes microscopically accessible due to an interplay of coherent and dissipative dynamics. To this end, in section 2 we work out in detail the underlying model and discuss its improvements in comparison to [38,39]. Based on the corresponding Lindblad master equation, we derive the resulting equations of motion of expectation values of the relevant system operators. In section 3 we determine the realistic model parameters in relation to current experiments. We then proceed in section 4 to analyse the steady-state properties of the system and identify the two limiting cases: a photon Bose-Einstein condensate and a laser-like regime. The latter one is novel and accessible precisely due to the inclusion of the coherent dynamics. For the former case, in section 5 we determine the dye-mediated dimensionless photon-photon interaction strength from realistic experimental data and, in particular, how it depends on the respective system parameters. Section 6 presents our concluding remarks.

Model
Let us now introduce a physical system which encompasses both laser and photon BEC as the possible limiting cases. This is going to be done in close correspondence with the actual experimental setups of photon BEC experiments. We consider N identical noninteracting two-level systems (TLS) inside an optical cavity. The transition between the two levels has the frequency ∆ and it is nearly resonant with M modes of the cavity. The dipole coupling between the TLS and the cavity modes has the strength g and it is assumed to be sufficiently weak so that the rotating wave approximation (RWA) holds. In the photon BEC experiments, the TLS were actually dye molecules dissolved in a solvent. The dye molecules have very broad rovibrational absorption and emission spectra, which can be modelled as an on-site phonon coupled to its own thermal bath [38,39]. In addition, due to frequent collisions with the solvent particles the dye molecules experience rapid dephasing. Hence, we take that each of the TLS is coupled to its own reservoir of R ≫ 1 harmonic oscillators. This can be thought of as a compound reservoir consisting of a phonon and its bath. The reservoirs are supposed to be independent and of identical properties. The collisional dephasing rate of each TLS is denoted by γ φ . We also assume that the TLS are incoherently pumped to the excited state with the rate γ ↑ and decay to the ground state with the rate γ ↓ via spontaneous emission of photons outside of the cavity. The decay rate of all cavity modes is abbreviated by κ. A conceptually similar system has previously been treated by Kirton and Keeling [38,39] using a mixture of the master equation and the Schwinger-Keldysh formalisms, but without accounting for the dephasing quantitatively. Our approach, instead, is based entirely on the master equation formalism and we improve several aspects of their model. Later on we underline those specific points and our enhancements that enable us to have access to a completely different regime of physical parameters.
In reality, the coupling between TLS and some cavity mode will also depend on the spatial mode function. A tractable model that incorporates the spatial dynamics was devised by Keeling and Kirton [40]. It has led to the successful understanding of the recent experiments [6,17]. However, the spatial dynamics introduces yet another level of complexity to the theoretical description. It could be implemented in our approach as well, but that would make the numerical calculations an order of magnitude more challenging. Thus, in the present work we make two additional simplifying assumptions: (i) all TLS are at exactly the same position and (ii) all cavity modes have the same intensity at the position of the TLS. This means that all TLS can be considered to evolve in an equivalent manner. Later on we will indicate how these assumptions may influence some of our results.

Master equation
Due to the above mentioned assumptions we consider the system Hamiltonian ( = 1) where ω m denote the cavity-mode frequencies and a m (a † m ) the bosonic annihilation (creation) operators of the cavity modes. The Hamiltonian H (j) ↑↓,R describes the j-th TLS and its reservoir, with σ ± j and σ z j being its Pauli spin operators. Bosonic annihilation (creation) operators and frequencies of the reservoir oscillators are b j,r (b † j,r ) and w r , respectively, while λ r are the appropriate interaction strengths. Since the experimental spectra of the dye molecules are very broad, we are led to assume that the TLS-reservoir coupling is strong. In order to treat it non-perturbatively to all orders, we perform the polaron transformationH = UHU † with and findH up to constant terms, where are the polaron displacement operators of the j-th TLS. In this way,Ṽ ↑↓,C captures the coupling of the TLS and the cavity modes which is dressed by the reservoir oscillators.
In order to proceed further, we assume that the oscillators in the polaron frame represent a bath in a thermal state at temperature T where Z β stands for the canonical partition function and β = 1/(k B T ) [66]. We consider such initial conditions that the subsystem TLS-cavity is uncorrelated with the bath in the polaron frame, i.e., ρ total = ρ ↑↓,C ⊗ ρ β . Since the coupling strength g is supposed to be weak, the bath influence can be incorporated by means of a master equation, i.e., by treatingṼ ↑↓,C as a perturbation up to the second order [67,68]. The first order contributes to the coherent unitary evolution through the thermal-averaged term where we introduced the notation X β ≡ Tr[Xρ β ] for a bath expectation value. Using the result for a harmonic oscillator of frequency ω in a thermal state, we find for the bath expectation value of the polaron displacement operators (4) Hence, one can naturally introduce a bath-dressed TLS-cavity coupling strength g β = g D ± j β . Obviously, due to (8) we have 0 < g β /g < 1, so that the influence of the bath in the first order is to effectively reduce the TLS-cavity interaction At this point we note that the previous first-order term was omitted by Kirton and Keeling [38,39], based on the implicit assumption that it is irrelevant due to the rapid collisional dephasing [69]. As we will demonstrate below, its influence deep in the photon BEC regime turns out to be negligible, so this regime can be described satisfactorily even if it is not taken into account. However, in the opposite laser-like regime such a term does play a major role, even in the presence of a fast sub-picosecond dephasing. Anyhow, on formal grounds, it should be a part of the proper treatment. We continue by applying the Born-Markov approximation as well as RWA, by tracing out the bath degrees of freedom and by taking into account the cavity losses along with the pumping and the decay of the TLS, similarly as [38,39]. As already mentioned, we additionally account for the dephasing of the individual TLS. Incoherent pumping can be formally described as coupling each TLS to a bath of inverted harmonic oscillators [70]. With this we find that the reduced density matrix ρ ↑↓,C of the TLS-cavity subsystem obeys the following master equatioṅ where L[X]ρ = {X † X, ρ} − 2XρX † and we have moved into the frame rotating with the frequency ∆, so that δ m = ω m − ∆ stands for the detuning of the cavity mode from the TLS transition. The thermal fluctuations ofṼ ↑↓,C give rise in the second order of perturbation theory to the incoherent transitions described by the dissipative Lindblad terms contained in the last double-sum of (10). The terms proportional to γ + m correspond to the absorption of the cavity photons by the TLS, while those with the prefactor γ − m represent the stimulated emission into the cavity modes. The previous approach should be satisfactory wheneverṼ ↑↓,C has small fluctuations around its thermal average and when the characteristic time scale, in which the bath modes undergo a displacement in order to adjust themselves to the instantaneous state of the TLS-cavity subsystem, is very short in comparison with the time scale of the subsystem relaxation. Note that the additional Lamb shifts due to the presence of the bath have been neglected as in [38,39]. Due to the dynamical influence of the bath, the corresponding rates γ ± m = γ(±δ m ) turn out to be frequency dependent and are obtained along the lines of [38,39,71] as with being the retarded connected correlation function of the bath displacement operators. One can notice that the pumping and the decay of the TLS yield an exponentially decaying factor in (11), i.e., they introduce an additional level broadening [71]. The time evolution of D − j (t) is generated by the free Hamiltonian H Note that in the long-time limit the many oscillatory terms from the above sum simply add up to zero, such that, recalling (8), one finds lim , the two displacement operators of very distant moments in time become uncorrelated. In the Kirton-Keeling's approach [38,39], the definition of the quantity (11) was actually without the last term of (12), i.e., C β (t) had a finite longtime limit. On formal grounds, if both γ ↑ and γ ↓ are zero, that leads to a divergence of the absorption and emission rates of resonant light, i.e., γ(δ = 0) → ∞. We trace this shortcoming back to the very absence of the first-order coherent term (9) from their treatment. Thus, the two improvements we have made to their approach come in pair.
The full master equation (10) is notoriously difficult to solve. However, its structure already reveals some general features of the system dynamics. Namely, one can clearly distinguish the coherent and the dissipative influence of the oscillator bath. The former one comes from the TLS-cavity coupling of the reduced strength g β -in a typical laserlike fashion, while the latter one is realised through the terms containing γ ± m , which were shown to lead to thermalisation of light and emergence of photon BEC [38,39]. In the following we will demonstrate that precisely their interplay determines these two limiting stationary behaviours, i.e., a photon BEC or a laser.

Equations-of-motion approach
In order to be able to perform a quantitative analysis, we proceed to obtain the equations of motion for the mean values of the system observables X ≡ Tr[Xρ ↑↓,C ], e.g., the populations of the cavity modes, the population inversion of the TLS etc. from the master equation (10). Since this procedure yields an infinite hierarchy of coupled equations, we use the cumulant expansion method [72][73][74][75] to truncate the hierarchy at the second level, i.e., we will keep the cumulants up to the second order only. If one wants to calculate higher-order correlation functions, a higher level of truncation will be necessary. However, due to the presence of coherent terms in the master equation, the situation becomes considerably more involved than in [38,39], even at this second-order truncation level.
We note that the system possesses a U(1) gauge symmetry: In a single experimental run a coherent field with a particular spontaneously chosen phase φ can build up. However, since the density matrix describes an average over many such realisations, the following equalities hold a m = a † m = σ − j = σ + j = 0 and similarly for all other gauge non-invariant operators [73]. In particular, it means that a m σ + where the index c denotes connected correlation functions. For instance, one has Due to the two simplifying assumptions mentioned in the beginning, we assume that all the TLS are mutually equivalent such that, e.g., σ z i = σ z 1 and σ + i σ − j = σ + 1 σ − 2 for all i, j = 1, . . . , N and i = j. The resulting equations of motion for the cavity mode whereas for the TLS population inversion σ z 1 we obtain d dt These equations are exactly like those of Kirton and Keeling [38,39], apart from the coherent terms proportional to g β which introduce the additional coupling to the mixedtype terms a m σ + 1 = a † m σ − 1 * . They measure the correlation between TLS and cavity photons and evolve according to where now the additional quantities a † k a m and σ + 1 σ − 2 appear. The former represent the correlations between different cavity modes, whose evolution is governed by and the latter the correlations between dipoles of different TLS following from It is important to notice that the quantities a m σ + 1 , a † k a m and σ + 1 σ − 2 can reach nonzero stationary values precisely due to the coherent part of the evolution, which we have introduced in addition to [38,39].
At this point, one additional specialisation is in order. Namely, based on the photon BEC experiments, we consider the photon modes as being transverse, arising from a twodimensional effective harmonic potential [5]. We thus consider regularly spaced cavity levels ω ℓ = ω 1 + (ℓ − 1)Ω with ℓ = 1, . . . , L, such that the energy level ω ℓ has the degeneracy d ℓ = 2ℓ, where the factor 2 comes from the two independent polarisations of light. The lowest frequency ω 1 represents the cavity cutoff. Since degenerate cavity modes evolve in the same manner, we have for any arbitrary function f , where from each level ω ℓ we have chosen a representative mode described by a ℓ and a † ℓ .

Bath model
In the following we analyse the stationary solutions of the equations (15)- (19) in the two regimes: a photon BEC and a laser-like regime. To this end, we specialise the model by choosing the bath spectral density, defined by J(w) = R r=1 λ 2 r δ(w − w r ), to be super-ohmic with an exponential cut-off [76,77] where η represents a dimensionless parameter measuring the coupling strength of the TLS to the bath and w c is the cutoff frequency of the bath. This allows us to obtain for (8) and (13) the closed-form expressions where ψ(z) = Γ ′ (z)/Γ(z) denotes a logarithmic derivative of the gamma function.
In relation to the experiments, both parameters η and w c characterise the impact of the used solvent on the spectral properties of the dye molecules. Below we discuss in detail that they change the absorption and the emission spectra γ(±δ) drastically. Furthermore, they turn out to have an impact upon the resulting bath-dressed coupling strength g β .

Determination of realistic model parameters
In order to apply our theory to the current experimental setups, we have to fix the model parameters in the experimentally accessible regimes. One of the key ingredients for the thermalisation of photons are the spectral properties of the dye molecules [17,50].
where the spectral (rovibrational) temperature T of the dye molecules is T = 300 K [17]. Comparing the rate equation from the Supplemental Material of [17] with our equation (15), we can interpret the rates γ + m = γ(δ m ) and γ − m = γ(−δ m ) as the absorption and the emission rates B 12 (ω m ) and B 21 (ω m ), respectively. Therefore, we fit the expression γ(ω − ∆) from (11), using equations (22a) and (22b) as well, to the experimentally measured absorption spectrum B 12 (ω) of the used Rhodamine 6G dye dissolved in ethylene glycol [81], by taking into account the absolute value B 12 (ω = 3300 THz) = 1.3 kHz [17]. The fitting allows then to fix the parameters of our bath model, namely the dimensionless coupling strength η of the TLS and the bath, the cutoff frequency of the bath w c , as well as the TLS-cavity coupling strength g and the TLS transition frequency ∆, which physically corresponds to the zero-phononline frequency of the dye. The obtained values are listed in table 1. Thus, the fitted value for η = 0.6 leads to a rather small bath-dressed TLS-cavity coupling strength g β = 8.3 × 10 −3 g, as expected, since this corresponds to the BEC regime. Figure 1 shows the resulting fit for γ(ω − ∆) and the experimentally provided data. Note, that we do not fit the absorption curve for ω > 3700 THz, since most of the relevant cavity modes are not influenced by the departure from the actual spectrum, as is indicated by the vertical dashed lines. If one wants to also fit this higher-frequency region, the inclusion of another bath would be necessary since the experimental spectrum displays the presence of another peak. A corresponding physical motivation in terms of another dye-molecule active phonon coupled to a thermal environment was discussed in [39]. With the fitted values our theory provides the emission rate curve γ(−ω + ∆) as well, which is shown in the same figure. The curves γ(ω − ∆) and γ(−ω + ∆) cross at the frequency ∆. We checked that the Kennard-Stepanov relation (23) is valid in the BEC regime for the relevant range of the cavity mode frequencies. For comparison, in the laser regime where η is small, the absorption and emission curves become squeezed towards the zero-phonon-line frequency ∆, see the dotted curves in figure 1. However, the Kennard-Stepanov relation (23) in this regime is no longer granted for frequencies highly detuned from ∆. The corresponding experimental number of dye molecules is taken to be N = 10 9 , based on the extensive discussion in [64]. The loss rate γ ↓ can Table 1. Parameters of the model adjusted to current experimental setup of photon BEC [5]. The first four parameters are obtained from the fit to the absorption spectrum of the dye. The remaining parameters are taken from [17,64,81]. be fixed to 0.25 GHz [81]. Furthermore, the pumping of the TLS γ ↑ is considered as a control parameter which is tuned to cross the boundary of the phase transition.
Knowing that dye molecules experience at least 10 12 collisions per second with solvent molecules [82], we will consider here γ φ = 0.1 THz as a referent value. Since there is an uncertainty about the exact order of magnitude of this rate, we will later analyse how its variation in a broad range of values influences our results.
In the experiment the cavity cutoff wavelength, corresponding to the mode of frequency ω 1 , can be tuned from 570 nm to 610 nm. The frequency separation of the cavity modes is 0.26 THz [5]. For the simulation we choose ω 1 to correspond to 585 nm, thus we get for the highest detuning δ 1 = ω 1 − ∆ = −260 THz. Our simulations cover the same spectral range as in [5], but due to computational complexity we choose a higher value Ω = 1.4 THz, if not stated otherwise. The photon loss rate κ is frequencydependent [17], so we take a mean value κ = 3.5 GHz. For the sake of clarity the used spectral range is marked in figure 1 by vertical dashed lines.

Two regimes: photon BEC and laser-like state
Having discussed the realistic values of the model parameters, in this section we take: N = 10 9 , w c = 20.5 THz, T = 300 K, g = 2.3 GHz, γ ↓ = 0.25 GHz, γ ↑ = 0.1 GHz, κ = 3.5 GHz and δ 1 = −260 THz. The results will be presented for two values of the dephasing rate, γ φ = 0.1 THz estimated from the literature [82] and a much larger one, γ φ = 10 THz, for the sake of comparison. Here we take Ω = 10 THz and consider L = 20 energy levels of the cavity. In the following we distinguish two limiting regimes: (i) η 1. In this case one has g β /g = D ± j β ≪ 1, meaning that the coherent contribution of the bath is highly suppressed and the evolution is dominated by the dissipative influence which leads to a thermalisation of light and an emergence of photon BEC; (ii) η ≪ 1. Here one finds g β /g ≈ 1, so that the bath has a pronounced coherent influence. In addition, the rates γ ± m of the highly detuned cavity modes acquire orders of magnitude smaller values in comparison with the previous regime, so that the dissipative influence becomes overwhelmed, but is still relevant. Hence, we expect that the non-equilibrium stationary state is then highly coherent and laser-like.
For the regime (i) we choose η = 1.5, yielding g β /g = 6.2 × 10 −6 . The resulting steady state solution yields the distribution of occupations of the cavity modes n ℓ for ℓ = 1, . . . , 20 as shown in figure 2(a). It is noticeable that the results (almost) do not depend on the value of dephasing rate. This is predictable since the coherent evolution is anyhow largely suppressed in this regime. The lowest energy level is macroscopically occupied with about 2.84×10 7 photons. The straight line corresponds to a Bose-Einstein distribution log(1 + 1/n ℓ ) = β(δ ℓ − µ), where µ denotes the chemical potential. Such a state was already analysed in detail in [38,39] and it was shown to correspond to a photon BEC. Our approach enables us to additionally characterise the stationary states by their photonic correlations. Namely, we have access to the quantities which provide a measure of correlations between representative cavity modes related to different energy levels. Their values belong to the interval [0, 1], where values close to 1 (0) correspond to a high (low) degree of correlation. In case (i) we find c ℓ,ℓ ′ < 4 × 10 −7 , i.e., the photon BEC state has almost no correlation between the modes of different frequencies. This is expected since the correlations build up through the coherent evolution which is highly ineffective in this regime.
In the opposite regime (ii), we take η = 0.05, which gives g β /g = 0.67. The distribution of stationary populations of the representative cavity modes is presented in figure 2(b) for two values of γ φ , namely 0.1 THz and 10 THz. In the former case, the cavity level ℓ = 9 acquires macroscopic occupation of almost 1.90 × 10 7 photons, but the distribution is quite distinct from the Bose-Einstein one. In this case we find c ℓ,ℓ ′ > 0.996, which demonstrates that the stationary state contains a quite high degree of correlations among the cavity modes of different energies. This is expected since the coherent influence of the bath is very pronounced. The stationary state is laserlike with some dissipative bath influence. For the larger value of the dephasing rate, the level ℓ = 10 becomes macroscopically occupied with 1.76 × 10 7 photons, while c ℓ,ℓ ′ > 0.986. Interestingly, in the considered parameter regime the results for γ φ = 0 would be almost indistinguishable from those for γ φ = 0.1 THz. This means that there is a certain dephasing threshold below which the coherent system dynamics is robust to the dephasing. Moreover, if one considered the full spatial structure of the cavity modes as in [40], the aforementioned correlations would decrease due to only partial overlap among different modes. However, this would not alter the present conclusions. We note that similar states supporting macroscopic occupations of optical modes of higher energies were also observed in [38,39] and, indeed, we can also reproduce such behaviour in the regime (i). However, the steady-state we have just presented features near-unity correlations of light, which represents a crucial difference and indicates that it is of an entirely different nature. Moreover, for different values of the cavity decay rate even the lowest level can acquire macroscopic occupation, while the populations of other cavity levels strongly depart from the Bose-Einstein distribution. Hence, the stationary states in the two regimes (i) and (ii) differ completely regarding the correlations and the distribution of photons among the cavity levels, as a consequence of different influences of the bath.

Properties of photon BEC
In the following, we focus on interesting properties of the photon BEC. At first, we determine from our model the equation of state. This allows then to extract the dimensionless effective photon-photon interaction strength in the photon BEC regime and study its dependence on various model parameters, which could be tuned experimentally. We will adopt the terminology in accordance with the experiments and, for instance, instead of TLS refer to dye molecules.

Equation of state
In this section we apply our theory to experimentally realistic values and determine at first the steady state of the equations of motion (15)- (19) for different pumping rates γ ↑ and evaluate the dependence of the chemical potential µ on the total photon number n tot , i.e., we obtain the equation of state µ(n tot ). In the following we present and discuss this procedure in detail by using the specific parameters from table 1, if not stated otherwise.
In figure 3(a) we show the occupation of different energy levels of the cavity in the BEC regime for the increasing pump rate γ ↑ with all other parameters being fixed. The onset of the condensation starts at the critical value γ cr ↑ . A zoom around γ cr ↑ is shown in figure 3(b). Clearly, the occupation of the lowest level d 1 n 1 shows a sudden increase at the critical point and becomes macroscopic afterwards. Further increase of γ ↑ mainly populates the lowest level, while the population of higher energy levels does not change significantly after the transition. At the critical value of the pumping parameter γ cr ↑ the total number of photons n tot = L ℓ=1 d ℓ n ℓ amounts to n tot (γ cr ↑ ) ≈ 2800, which is quite close to the expected critical photon number for the condensation onset in the case of non-interacting interacting bosons in two dimensions [15,16] and at T = 300 K where the existence of two independent polarisations of light has been taken into account. Note that the rise of the total photon number n tot around the transition point becomes smoothened when the number of cavity levels L in the considered frequency range is increased, as is seen in figure 3(b). This can be explained with the significant contribution of degeneracies d ℓ of the energy levels with high ℓ to n tot . The mode occupation n ℓ for a fixed γ ↑ is shown in a logarithmic representation in figure 4(a), along the triple-dashed vertical line of figure 3(a). The linear behaviour indicates that the modes are distributed according to Bose-Einstein statistics Fitting a linear function to log (1 + 1/n ℓ ) yields the inverse temperature β and the chemical potential µ due to equation (26). For a non-interacting BEC condensate, i.e., g β = 0, the temperature obtained by fitting to a thermal cloud coincides with the spectral temperature of the dye and the chemical potential is locked to the value δ 1 above the threshold [38]. However, here g β is non-zero but small, which induces additional small corrections of the occupation. As a consequence, the effective thermalisation temperature of the thermal cloud differs generically from the spectral temperature of the dye. This behaviour depends on all system parameters, and especially a high pumping rate γ ↑ can significantly change the temperature. For a choice of parameters outside of a certain region the non-equilibrium properties of the model are dominating the tendency to thermalise and the distribution is, in general, not thermal any more. This happens even in the case g β = 0 [39]. Therefore, we restrict our further analysis only to the cases where a thermal cloud does exist. We have observed that a few of the lowest modes do not follow the linear dependency in figure 4(a), thus they are not considered in our fitting procedure. Additionally, we drop some of the highest modes as well, as they hold relatively small occupation and can, therefore, have large relative numerical error. The resulting value of µ from figure 4(a) is shown by a star symbol in figure 4(b). We repeat the procedure for different values of the pumping rate γ ↑ and obtain a linear equation of state µ = µ(n tot ), as presented in figure 4(b).

Photon-photon interaction strength
The slope ∂µ/∂n tot of the equation of state µ = µ(n tot ) in figure 4(b) is a consequence of the dye-mediated effective photon-photon interaction (see Appendix A), which is measured by the dimensionless interaction strength Thus, we infer from figure 4(b) the resulting interaction strengthg = 5.2 × 10 −7 , which agrees quite well with the rangeg ∼ 10 −8 − 10 −7 given in [63]. Note that, as mentioned in the introduction, the thermo-optical effect dominates the photon-photon interaction. Therefore, this value cannot be directly compared with the measured values [5,62,65], since we only model one of the respective contributions.  (see table 1), apart from L = 101 which is used in (b)-(h). As we see, increase of N or decrease of T significantly increasesg, whereas κ and γ φ affect its value only slightly in the vicinity of experimentally realistic values. Fit (red line) in (e) indicates the dependence ofg on g β in the form of an even polynomial. Now we investigate howg depends on the parameters of our model. First, we vary those that could be tuned in experimental setups: the number of cavity levels L in the experimentally relevant range of 140 THz above δ 1 , the number of dye molecules N, the spectral temperature of the dye T and the cavity decay rate κ. Next, it is instructive to think of g β as being an independent parameter. In this way, by tuning g β from zero to its experimental value g exp β , we are able to examine how the coherent terms of our model give rise to the effective photon-photon interaction strength. Finally, since the dephasing rate γ φ is only approximately known in the experiments, we also investigate its influence ong for various values of other parameters. The corresponding results are presented in figure 5. The number of equations increases quadratically with the number of levels L. We have chosen L = 101 due to computational constraints. Figure 5(a) shows that the interaction strengthg depends non-monotonously on L. In case of up to 200 levels, the interactiong increases nearly exponentially with the number of cavity levels. However, after reaching the maximum at L ≈ 200 the dimensionless interaction strengthg starts to decrease. The value ofg in the case of the experimental number of levels L = 501 is then comparable with the value for L = 101 levels, which we used in our simulations.
According to figure 5(b), a larger number of dye molecules N increasesg dramatically, since more dye molecules mediating an effective coupling between the photons are present. In addition, much larger photon-photon interaction can also be achieved by lowering the temperature, see figure 5(c). In contrast, figure 5(d) reveals that increasing the decay rate κ decreases only slightly the interaction strengthg. An intuitive explanation of the last two results is as follows. Increasing either the temperature or the photon decay rate reduces the total number of photons in the system, thus there are less photons available to modify the dye medium andg decreases correspondingly.
In figure 5(e) we investigate the dependency ofg on the coherent coupling g β , which we vary artificially from zero to the experimental value g ext β . Quite expectedly, in the limit g β → 0 the dimensionless interaction strengthg practically vanishes, the latter corresponding to the case of the Kirton-Keeling model [38,39]. The increase of g is non-linear and in the considered range an even polynomial in g β yields a good fit (red line). We observe that the quadratic term is almost negligible compared to quartic and higher-order terms, akin to the consideration of a photon-photon interaction corresponding to the box Feynman diagram analysed in the work [63]. In our framework, such a dependence could easily be understood via a simple perturbative expansion of the expectation values X = ∞ p=0 X (p) , where X (p) ∝ g p β , for an arbitrary system operator X in the equations (15)- (19), with the time derivatives set to zero. In such a way, higher perturbation orders can be systematically calculated from the lower ones. The zeroth order solutions for n m (0) and σ z 1 (0) are exactly those from the work of Kirton and Keeling [38,39]. It turns out that n m = n m (0) + n m (2) + n m (4) + . . ., so that we have an expansiong =g (2) +g (4) + . . . in even powers of g β .
Next, we analyse the effect of the dephasing ong and the sensitivity of the obtained dependence on other system parameters. Figures 5(f)-5(h) show that, quite generically, g is affected insignificantly when the dephasing rate γ φ is varied from zero to few THz. Such a result could be attributed to the presence of a large detuning in (17), of the order of 100 THz, which anyhow strongly suppresses the coherent evolution on its own. Further increase of γ φ may lead to the appearance of a resonance-like peak, followed by a polynomial decay for dephasing rates above several hundreds of THz. The latter is expected to happen when the dephasing rate becomes larger than |δ 1 |. The resonancelike peak becomes more pronounced as g β is increased or the temperature decreased, as is seen in figures 5(f) and 5(h). This observation indicates that the peak is of coherent origin. On the other hand, the resonance-like peak can completely disappear when the cavity cutoff frequency is shifted towards the zero-phonon line, which is demonstrated in figure 5(g).

Conclusions
Here we presented a model which can interpolate between two different kinds of states of light in a micro-cavity, namely between a nearly non-interacting photon BEC and a laser-like state. Our model is based on a master equation approach, with an interplay between coherent and dissipative dynamics. The dominance of the former or the latter leads either to a coherent lasing state or to an equilibrated BEC state, respectively. We demonstrated that in the BEC case the lowest cavity energy level is macroscopically occupied and cavity modes of different energies are almost uncorrelated, whereas in the lasing case some cavity level becomes macroscopically occupied and strongly correlated with the others. Afterwards, we showed how to fix the parameters of our theory in an experimentally realistic regime. We emphasised that the coherent part of the master equation is then overwhelmed by the dissipative effects, but still large enough to lead to an additional effective photon-photon interaction. As a consequence, the chemical potential depends linearly on the total number of photons, as is expected from a perturbative solution of a Gross-Pitaevskii equation for the condensate wave function. This dependency allowed us to determine the dimensionless interaction strengthg to be of the order of 10 −7 for experimentally realistic parameters. We also investigated the dependency ofg on different model parameters, which can feasibly be tuned in the photon BEC experiments. Our numerics showed that increasing the number of dye molecules or decreasing the spectral dye temperature can significantly increase theg value, whereas it is not much influenced by the cavity loss rate κ. However, this value cannot be directly connected to the current experimental values [5,62,65]. The reason is that in the experimental setups the dominating photon-photon interaction is of thermo-optical origin, whereas our theory has no spatial degrees of freedom and, thus, cannot capture such diffusive effects. Instead, in our case the effective photonphoton interaction could be compared with the dye-mediated photon-photon scattering. And indeed, our value is in the range of the expected estimate [63].
Another currently disputed feature of the photon condensate concerns its possible polarisation. Whereas no significant polarisation of the photon BEC was found in the original Bonn experiment [5], recent systematic measurements of the Stokes parameters in Utrecht [83] indicate that the polarisation of the photon BEC correlates with the polarisation of the pump pulse. These new experimental results together with the recent theoretical investigation in the BEC case [41] offer the prospect that the polarisation dependency could be investigated on the basis of an extension of our microscopic model during the whole crossover from the photon BEC to the laser-like phase.
With this ansatz the Gross-Pitaevskii equation (A.1) reduces to which determines both interaction corrections Ψ (1) (x) and µ (1) . In our context it is sufficient to calculate the latter one, which follows from the Fredholm alternative [84].
To this end we multiply equation (A.5) with Ψ (0) (x) and integrate over x, so we get due to (A.2) with the dimensionless interaction parameter [61]  Thus, for small interactionsg the chemical potential µ changes linearly with the photon number n tot . The slope ∂µ/∂n tot of the equation of state µ(n tot ) depends then via ∂µ/∂n tot =g Ω/(2π) linearly on the dimensionless interaction strengthg, which leads to relation (27).