Decondensation in non-equilibrium photonic condensates: when less is more

TL;DR

This paper explores the complex phase behavior of photons in a dye-filled microcavity, revealing how varying pump and thermalization can lead to diverse states like Bose-Einstein condensation, multimode condensation, lasing, and unexpected decondensation phenomena.

Contribution

It provides a detailed non-equilibrium phase diagram of photonic states, highlighting the novel decondensation of modes under certain conditions.

Findings

Identification of multiple phase transitions including condensation and lasing

Observation of mode decondensation in non-equilibrium conditions

Rich phase diagram illustrating complex photon behaviors

Abstract

We investigate the steady state of a system of photons in a pumped dye-filled microcavity. By varying pump and thermalization the system can be tuned between Bose-Einstein condensation, multimode condensation, and lasing. We present a rich non-equilibrium phase diagram which exhibits transitions between these phases, including decondensation of individual modes under conditions that would typically favor condensation.

Equations

$$\begin{array}{rl}\frac{d{n}_m}{dt}=-\kappa n_m+& \rho {\Gamma }_{\downarrow }^{(m)}{f}_m(n_m+1)\\ +& \rho {\Gamma }_{\uparrow }^{(m)}(f_m-1)n_m,\end{array}$$

$$\frac{\partial f(\mathbf{r})}{\partial t}=-{\Gamma }_{\downarrow }^{tot}(\mathbf{r})f(\mathbf{r})+{\Gamma }_{\uparrow }^{tot}(\mathbf{r})(1-f(\mathbf{r})),$$

$$\Gamma^{tot}_{k}(\mathbf{r})=\Gamma_{k}(\mathbf{r})+\sum_{m}|\psi_{m}(\mathbf{r})|^{2}\Gamma_{k}^{(m)}(n_{m}+\delta_{k\downarrow})\$$

$$\begin{array}{r}{f}_m^c=\frac{{\Gamma }_{\uparrow }^{(m)}+\frac{\kappa }{\rho }}{{\Gamma }_{\downarrow }^{(m)}+{\Gamma }_{\uparrow }^{(m)}}.\end{array}$$

$$f_m^{(s)}=\int d^2\mathbf{r}|{\psi }_m(\mathbf{r}){|}^2\frac{{\Gamma }_{\uparrow }^{tot}(\mathbf{r})}{{\Gamma }_{\uparrow }^{tot}(\mathbf{r})+{\Gamma }_{\downarrow }^{tot}(\mathbf{r})}$$

$$n_m=\frac{{\Gamma }_{\downarrow }^{(m)}{f}_m}{\frac{\kappa }{\rho }-{\Gamma }_{\downarrow }^{(m)}{f}_m-{\Gamma }_{\uparrow }^{(m)}(f_m-1)}.$$

$$f_m^{(c)}=\frac{{\Gamma }_{\uparrow }^{(m)}+\frac{\kappa }{\rho }}{{\Gamma }_{\uparrow }^{(m)}+{\Gamma }_{\downarrow }^{(m)}}$$

$$\begin{array}{rl}{\Gamma }_{\downarrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\downarrow },\\ {\Gamma }_{\uparrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\uparrow }.\end{array}$$

$$f_m^{(s)}=\int d^2\mathbf{r}|{\psi }_m(\mathbf{r}){|}^2\frac{{\Gamma }_{\uparrow }}{{\Gamma }_{\uparrow }+{\Gamma }_{\downarrow }}=\frac{{\Gamma }_{\uparrow }}{{\Gamma }_{\uparrow }+{\Gamma }_{\downarrow }}.$$

$$\frac{{\Gamma }_{\uparrow }}{{\Gamma }_{\downarrow }}=\frac{{\Gamma }_{\uparrow }^{(m)}+\frac{\kappa }{\rho }}{{\Gamma }_{\downarrow }^{(m)}-\frac{\kappa }{\rho }}.$$

$$\begin{array}{rl}{\Gamma }_{\uparrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\uparrow }(\mathbf{r})+|{\psi }_{[0,0]}(\mathbf{r}){|}^2{\Gamma }_{\uparrow }^{[0,0]}{n}_{[0,0]}\\ {\Gamma }_{\downarrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\downarrow }(\mathbf{r})+|{\psi }_{[0,0]}(\mathbf{r}){|}^2{\Gamma }_{\downarrow }^{[0,0]}{n}_{[0,0]}.\end{array}$$

$$f_m^{(s)}=\int d^2\mathbf{r}|{\psi }_m(\mathbf{r}){|}^2\frac{{\Gamma }_{\uparrow }(\mathbf{r})+|{\psi }_{[0,0]}(\mathbf{r}){|}^2{\Gamma }_{\uparrow }^{[0,0]}{n}_{[0,0]}}{{\Gamma }_{\updownarrow }(\mathbf{r})+|{\psi }_{[0,0]}(\mathbf{r}){|}^2{\Gamma }_{\updownarrow }^{[0,0]}{n}_{[0,0]}},$$

$$f_{[0,0]}^{(s)}=\frac{{\Gamma }_{\uparrow }^{[0,0]}+\frac{\kappa }{\rho }}{{\Gamma }_{\uparrow }^{[0,0]}+{\Gamma }_{\downarrow }^{[0,0]}},$$

$$\begin{array}{rl}{\left|{\psi }_{[0,0]}\right|}^2=& \frac{1}{\pi }\frac{1}{L^2}{e}^{\text{nicefrac}-(x^2+y^2{)}^2{L}^2}\\ {\left|{\psi }_{[0,1]}\right|}^2=& 2{\left|{\psi }_{[0,0]}\right|}^2\frac{x^2}{L^2}\\ {\left|{\psi }_{[0,2]}\right|}^2=& \frac{1}{2}{\left|{\psi }_{[0,0]}\right|}^2\left[4{\left(\frac{x}{L}\right)}^4-4{\left(\frac{x}{L}\right)}^2+1\right]\\ {\left|{\psi }_{[0,3]}\right|}^2=& \frac{1}{3}{\left|{\psi }_{[0,0]}\right|}^2\left[4{\left(\frac{x}{L}\right)}^6-12{\left(\frac{x}{L}\right)}^4+9{\left(\frac{x}{L}\right)}^2\right]\end{array}$$

$$\begin{array}{rl}{f}_{[0,0]}=& 2I_0\\ f_{[0,1]}=& 2I_2\\ f_{[0,2]}=& \frac{3}{2}{I}_4-2I_2+I_0\\ f_{[0,3]}=& \frac{5}{6}{I}_6-3I_4+3I_2,\end{array}$$

$$I_k=\frac{{\Gamma }_{\uparrow }}{{\Gamma }_{\updownarrow }}\int dr\frac{1+\sigma e^{-r^2}}{1+\eta e^{-r^2}}{r}^{(k+1)}{e}^{-r^2},$$

$$\begin{array}{rl}\sigma =& \frac{{\Gamma }_{\uparrow }^{[0,0]}}{{\Gamma }_{\uparrow }}\frac{n_{[0,0]}}{\pi },\\ \eta =& \frac{{\Gamma }_{\updownarrow }^{[0,0]}}{{\Gamma }_{\updownarrow }}\frac{n_{[0,0]}}{\pi }.\end{array}$$

$$\begin{array}{rl}{I}_0=& \frac{{\Gamma }_{\uparrow }}{2\eta {\Gamma }_{\updownarrow }}\left(\sigma +\left(1-\frac{\sigma }{\eta }\right)\log (1+\eta )\right)\\ I_2=& \frac{{\Gamma }_{\uparrow }}{2\eta {\Gamma }_{\updownarrow }}(\sigma \left(1-\frac{{\pi }^2}{6\eta }\right)-{\mathrm{Li}}_2(-\eta )+\\ & \frac{\sigma }{\eta }\left(\frac{\log (1+\eta )}{2}\log \left(\frac{1+\eta }{{\eta }^2}\right)+{\mathrm{Li}}_2\left(\frac{1}{1+\eta }\right)\right))\\ I_4=& \frac{{\Gamma }_{\uparrow }}{\eta {\Gamma }_{\updownarrow }}\left(\sigma -\frac{1-\sigma }{\eta }{\mathrm{Li}}_3(-\eta )\right)\\ I_6=& \frac{3{\Gamma }_{\uparrow }}{\eta {\Gamma }_{\updownarrow }}\left(\sigma -\frac{1-\sigma }{\eta }{\mathrm{Li}}_4(-\eta )\right)\end{array}$$

$${\mathrm{Li}}_j(x)=\sum _{i=1}^{\infty }\frac{x^i}{i^j}.$$

$$\begin{array}{rl}{\Gamma }_{\uparrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\uparrow }(\mathbf{r})+\sum _{p\in S_{0,2}}|{\psi }_p(\mathbf{r}){|}^2{\Gamma }_{\uparrow }^{(p)}{n}_p\\ {\Gamma }_{\downarrow }^{tot}(\mathbf{r})& \simeq {\Gamma }_{\downarrow }(\mathbf{r})+\sum _{p\in S_{0,2}}|{\psi }_p(\mathbf{r}){|}^2{\Gamma }_{\downarrow }^{(p)}{n}_p,\end{array}$$

$$f_m^{(s)}=\int d^2\mathbf{r}|{\psi }_m(\mathbf{r}){|}^2\frac{{\Gamma }_{\uparrow }(\mathbf{r})+\sum |{\psi }_p(\mathbf{r}){|}^2{\Gamma }_{\uparrow }^p{n}_p}{{\Gamma }_{\updownarrow }(\mathbf{r})+\sum |{\psi }_p(\mathbf{r}){|}^2{\Gamma }_{\updownarrow }^p{n}_p},$$

Full text

Phase transitions in systems governed by quantum statistics at thermal equilibrium have been investigated intensely. Canonical examples of these transitions are the formation of Cooper pairs in the superconductivity transition Bardeen1957 , or the transition of a thermal cloud of bosonic atoms to a Bose-Einstein (BEC) condensate as temperature decreases Anderson1995 ; Davis1995 . The study of non-equilibrium phase transitions has spanned many disciplines from classical physics to social sciences Lubeck2004 , and explained many phenomena, such as trends in society Weidlich1991 , or traffic jams Kerner1997 . Only recently has attention turned to non-equilibrium phase transitions in quantum systems. For example the condensation of polaritons in semiconductor Kasprzak2006 ; Deng2006 ; Balili2007 or organic Plumhof2014 ; Daskalakis2014 solids has been observed outside thermal equilibrium. Driven-dissipative many-particle quantum systems showing intricate phase diagrams PhysRevA.91.033609 ; PhysRevX.7.011016 , and bistability phenomena giving rise to coexisting phases have also been observed Sarchi2008 ; PhysRevX.7.011012 .

Our system of study is a gas of photons confined in a dye-filled microcavity which, when the dye is pumped, can be made to thermalize and Bose-Einstein condense Klaers2010 ; Klaers2010a , as predicted by thermal-equilibrium theory. Thermalization results from absorption and re-emission of light from the cavity by the dye Schmitt2015 , which is limited by emission from the cavity Marelic2015 ; Keeling2016 . Thermal equilibrium is thus always imperfect and breaks down completely if the cavity is far detuned from the dye molecular resonance Kirton2013 ; Kirton2015 . The system then features multimode condensation Marelic2016 ; Keeling2016 as a clear signature of non-equilibrium behavior. It has been noted that multimode systems driven far from equilibrium can show multimode condensation vorberg2013 , and that the kinetics of two-mode laser systems can be made to show a sort of minimalist Bose-Einstein condensation redlich2016 ; leymann2017 .

Extrapolating from prior experiments on non-equilibrium phase transitions PhysRevX.7.011016 ; PhysRevX.7.011012 , one would expect that condensation is always favored by an increase in the pump rate Kasprzak2006 ; Klaers2010 and a decrease in thermalization rate. This is typically true for the overall fraction of photons in condensed modes. For individual modes however, we also find the opposite behavior in the regime far out of equilibrium. That is, modes with condensed photons can loose their macroscopic occupation as the pump rate is increased, which is similar to decreasing entropy as temperature increases, or negative heat capacities michaelian2007 . This gives rise to a highly complex dependence of the steady state on properties like pump rate, geometry, and the time-scale of thermalization. Based on a microscopic model Kirton2013 we predict the dependence of multimode condensation on such properties. Despite the non-linear nature of the system we find analytic laws characterizing condensation and decondensation that coincide accurately with numerically exact solutions.

In dye-filled microcavities only one longitudinal mode of the harmonic cavity is sufficiently close to resonance with the dye molecules; the dynamics can then be reduced to a two-dimensional model with transverse modes labeled by $m=[m_{x},m_{y}]$. Taking into account loss through the mirrors and spontaneous emission from dye molecules into modes that are not confined within the cavity results in the well-established equation of motion Keeling2016

[TABLE]

for the average occupation $n_{m}$ of mode $m$, where $\kappa$ is the cavity decay constant, $\Gamma_{\uparrow}^{(m)}$ and $\Gamma_{\downarrow}^{(m)}$ are respectively the rates of absorption from and emission into mode $m$, and $\rho$ is the areal density of molecules; $f_{m}$ is the fraction of excited molecules interacting with mode $m$, and it is given in terms of the fraction $f(\mathbf{r})$ of molecules at point ${\bf r}$ which are excited and the mode-profile $\psi_{m}(\mathbf{r})$ via the relation ${f_{m}=\int d^{2}{\bf r}f(\mathbf{r})|\psi_{m}(\mathbf{r})|^{2}}$. The dynamics of the excited-state population is governed by

[TABLE]

in terms of the rates of total absorption and emission, $\Gamma^{tot}_{\downarrow}(\mathbf{r})$ and $\Gamma^{tot}_{\uparrow}(\mathbf{r})$, which depend on the mode occupations via

[TABLE]

with $k=\uparrow,\downarrow$ and $\delta_{\uparrow\downarrow}=1-\delta_{\downarrow\downarrow}=0$, where $\Gamma_{\uparrow/\downarrow}(\mathbf{r})$ are the pump rate of molecules by the laser, and the decay of molecules not captured by emission into the modeled cavity modes.

As the thermalization process occurs when excitations are exchanged between cavity modes and dye molecules, we compare the rate of absorption $\rho\Gamma_{\uparrow}^{(m)}$ to the rate of loss $\kappa$ and define the thermalization coefficient $\gamma=\rho\Gamma_{\uparrow}^{[0,0]}/\kappa$ 111The thermalization coefficient $\gamma$ is directly related to the detuning, but since this relation depends on the absorption spectrum of the molecules there is no simple analytical form for it. Nevertheless, since $\gamma$ is the more informative quantity, we will mostly refer to thermalization in the following.. We choose parameter values appropriate to real experiments, but specify all values in units of cavity decay $\kappa$ (which for typical experiments is of the order of $10^{9}/$s) and harmonic oscillator length $L$ (the mean spatial extent of the lowest cavity mode). The shape of absorption and emission profiles $\Gamma_{\uparrow/\downarrow}^{(m)}$ for the individual modes are extracted from experimental data nyman_robert_andrew_2017_569817 and their peaks are set to $1.2\times 10^{-9}\,\kappa$ (see appendix B). We consider a slightly anisotropic cavity with mode spacings $\omega_{x}=\omega_{y}/1.01=3\times 10^{4}\,\kappa$, and a detuning between the molecular resonance frequency and the lowest cavity eigenfrequency ranging from $-3.5\times 10^{5}\kappa$ to $-1.89\times 10^{5}\kappa$. The pump of the molecules of areal density $\rho=10^{12}/L^{2}$ has a Gaussian profile with width $20L$, and the decay rate $\Gamma_{{}_{\downarrow}}$ of excited molecules is set to $\Gamma_{{}_{\downarrow}}=\kappa\,/\,4$. The anisotropy of the cavity is chosen in order to avoid degeneracies, but it is sufficiently small so that mode pairs $[m_{x},m_{y}]$ and $[m_{y},m_{x}]$ behave almost identically, and their condensation thresholds are hardly distinguishable. We will therefore only discuss modes $[m_{x},m_{y}]$ with $m_{x}\leq m_{y}$. In order to arrive at a finite-dimensional problem, we consider cavity modes with $m_{x}+m_{y}\leq 6$ only.

Fig. 1 (top panel) depicts the stationary solutions for mode occupations as functions of pump rate, and one can see step-like increases and decreases of the populations, i.e. condensation and decondensation of individual modes at specific values of pump rate. The BEC phase is defined by condensation in the lowest cavity mode only, whereas a multimode condensate contains additional condensed modes. Any phase with one or more condensed modes, but an un-condensed ground mode is considered a laser 222The distinction we make here is that in Bose-Einstein and multimode condensate phases both absorption and emission play important roles, whereas (stimulated) emission dominates the behavior in the laser phase. and any phase without any condensed modes will be called un-condensed. With increasing pump rate, Fig. 1 (top panel) thus features the transition from an un-condensed phase to BEC, followed by a transition from BEC to a multimode condensate, and three transitions between different multimode condensates.

Despite the system’s complex, non-linear behavior, we can develop an understanding of (de)-condensation in terms of the decomposition of the total pump rate into individual contributions given in Eq. (3). For low populations $n_{m}$ of all modes, the total rates of absorption/emission can be well approximated by direct pumping and loss, i.e. Eq. (3) reduces to $\Gamma^{tot}_{\uparrow/\downarrow}(\mathbf{r})\simeq\Gamma_{\uparrow/\downarrow}(\mathbf{r})$, so that the pumping of dye molecules (Eq. (2)) is proportional to the external pumping. In the case of strongly occupied modes, on the other hand, the approximation $\Gamma^{tot}_{\uparrow/\downarrow}(\mathbf{r})\simeq\sum_{m}|\psi_{m}(\mathbf{r})|^{2}\,\Gamma^{m}_{\uparrow/\downarrow}n_{m}$ in terms of the rates $\Gamma^{m}_{\uparrow/\downarrow}n_{m}$ of stimulated absorption and emission holds, and the proportion of excited dye molecules is largely independent of external pumping but depends mostly on the interaction with the condensed mode or modes. Molecules in spatial domains where $|\psi_{m}(\mathbf{r})|^{2}\,\Gamma^{m}_{\uparrow/\downarrow}n_{m}\gg\Gamma_{\uparrow/\downarrow}(\mathbf{r})$ are thus clamped to mode $m$.

As pumping is increased the excitation of clamped molecules decreases, but the excitation of un-clamped molecules grows. In particular, un-condensed modes thus experience a growing reservoir of excitations that can help them to gain enough population to condense. As a mode condenses it also starts to clamp dye molecules and the number of molecules clamped to this mode grows with increasing population $n_{m}$. This mode thus enters a competition for access to excitations with the other modes, and potentially reduces the access for other modes, which can result in their decondensation. This effect can be explicitly seen in the fraction of excited molecules $f_{m}$ accessible to mode $m$ as shown in Fig. 1 (center and bottom panel). The center panel magnifies $f_{[0,1]}$ and shows that $f_{[0,1]}$ decreases with increasing pump rate after the condensation of mode $[0,2]$, despite the general (and expected) trend that $f_{m}$ grows with increasing pump rate. Quite strikingly, the decrease in accessible excited dye molecules required for decondensation is rather minute and visible only in the magnified panel. Nevertheless, it has an extremely significant impact and can result in decondensation of a mode.

With this qualitative understanding of decondensation at hand, we can now proceed to a quantitative prediction of condensation thresholds; a more detailed analysis and the explicit analytic solutions can be found in appendix A. The mode population $n_{m}$ for stationary solutions to Eq. (1) diverges if $f_{m}$ approaches the critical value

[TABLE]

Neglecting contributions from uncondensed modes to the total pump rate (Eq. (3)), one obtains the stationary solutions

[TABLE]

of Eq. (2) for all modes as functions of the mode populations $n_{m}$ of the condensed modes. Gain clamping implies the condition $f_{m}^{s}\lesssim f_{m}^{c}$ for the condensed modes, which, in turn determines $n_{m}$ for all condensed modes. This permits identification of pump rates (or other system parameters) that achieve the condensation condition for an originally uncondensed mode. The thresholds for Bose-Einstein condensation, multimode condensation and the decondensation of mode $[0,1]$ while mode $[0,2]$ is condensed are depicted as white lines in the phase diagram depicted in Fig. 2.

The accuracy of these estimates can be verified by comparisons to numerically exact solutions. The steady-state solution of Eq. (1), together with Eq. (5), can be found with algebraic root finding routines, and we verified explicitly that the obtained solutions coincide with the solutions obtained by propagation until a stationary state is reached. Defining phase boundaries can be done unambiguously in the thermodynamic limit, but it poses an intricate problem in systems of finite size. Since numerically exact solutions do not result in diverging mode populations, we employ large increases in the population of a single mode under small changes of driving conditions (as depicted in Fig. 1 (top panel)) as indicators of condensation threshold Kasprzak2006 ; Marelic2016 . Quantitatively, we define the threshold of condensation (decondensation) as an increase (decrease) in population greater than 3 orders of magnitude over an increase in pump rate of 10%, but due to the sharp nature of all observed transitions, the identification of threshold is largely independent of the explicitly chosen numerical values, see appendix C.

As depicted in Fig. 2 the system at lower pump rates is in the un-condensed phase, with no condensed modes. At large thermalization coefficient and pump rate the system is in the BEC phase. As the thermalization coefficient is decreased from here, the system passes through various multimode phases, until the ground mode is no longer condensed. Under these conditions the system behaves like a laser, where condensation is a consequence of stimulated emission, rather than a condensate, for which absorption and emission play equally important roles. The trend that strong pumping and weak thermalization favors condensation is observed for the onset of Bose-Einstein and multi-mode condensation, but in the regime of multi-mode condensation the behavior becomes more complicated and decondensation can be induced through an increase in pump rate $\Gamma_{\uparrow}$ or a decrease in thermalization coefficient $\gamma$. Given the extremely sensitive dependence of the system’s phase on the excitation of dye molecules, there are several instances where small changes in pumping or thermalization result in a substantial redistribution of populations between the modes, and multiple triple- and quadruple-points.

Eventually, a critical assessment of experimental relevance is in order. Current experiments are performed for $\gamma$ ranging between $0.2$ and $5$ and pump rates reach rates ten times the threshold value Marelic2016 . Exploring the upper half of the phase diagram presented in Fig. 2, requires substantially stronger pumping than currently realized. Since pulsed lasers can achieve peak pump rates three orders of magnitude higher than CW lasers, with pulse durations longer than the time taken to reach steady state, this seems a perfectly viable option. The model we have used is applicable to a wide range of systems which need only satisfy a few criteria: an optical environment with a well-defined ground state, a fluorescent gain medium, and re-scattering of light faster than loss resulting in thermalization. We thus expect that multimode condensation and decondensation both inside and outside the lasing regime will be observable in plasmonic lattices coupled to dyes, which have recently shown condensation hakala2017 , semiconductors in photonic crystal resonators deLeeuw16 , and also more conventional laser systems. Exciton-polariton condensates (both semi-conductor and organic), however, differ from photon condensates in that the mixed light-matter excitations relax directly through their matter component, whereas in photon and plasmon condensates equilibration occurs via exchanges between weakly-coupled light and molecular excitations. Our predictions therefore do not translate directly, but it would be intriguing to identify similar mechanisms in these condensates.

The rich interplay between the different modes observed in Fig. 2 also offers great opportunities for the creation of tailored states of light hisch2013 , since the mode structure can easily be influenced through the shape of the cavity flatten2015 , lattice hakala2017 or crystal structure deLeeuw16 , and the spatial pump profile can be varied. Since fluctuations are most relevant near phase transitions, we expect this rich phase diagram to be a fruitful tool in the search for unusual quantum correlations.

We acknowledge fruitful discussions with David Newman, Peter Kirton and Jonathan Keeling, and financial support from the UK EPSRC (via fellowship EP/J017027/1 and the Controlled Quantum Dynamics CDT EP/L016524/1), and the ERC (via ODYCQUENT grant). The data underlying this article are available at nyman_robert_andrew_2017_569817 .

Appendices

In appendix A we provide a detailed, analytic derivation of phase boundaries resulting from gain clamping. Appendix B contains a discussion of the absorption and emission spectra used in the simulation, followed by appendix C with a justification for the choice of threshold condition. In appendix D, we show how a variation of pump profile induces multimode condensation and decondensation, similarly to the variation of thermalization that is discussed in the main paper.

Appendix A Analytic Phase Boundaries

As sketched in the main paper, phase boundaries can be estimated analytically with high accuracy based on the clamping mechanism. Here we discuss the derivation of the boundaries depicted in Fig. 2 in more detail. In the presently considered case of spatially homogeneous pumping, both pump profile $\Gamma_{\uparrow}({\bf r})$ and molecular decay $\Gamma_{\downarrow}({\bf r})$ are constant and will be denoted by $\Gamma_{\uparrow}$ and $\Gamma_{\downarrow}$ in the following. We will also assume a rotationally invariant cavity with degenerate modes.

The population of mode $m$ in the stationary solution of Eq. (1) reads

[TABLE]

It diverges if $f_{m}$ reaches the critical density

[TABLE]

of excited molecules. Approximating Eq. (3) according to a given phase transition, and equating $f_{m}^{(c)}$ with Eq. (5), thus allows us to estimate phase boundaries.

A.1 Boundary between the un-condensed phase and single-mode condensation

The decay of the molecular excited states is dominated by spontaneous emission into free space at a rate $\Gamma_{\downarrow}$, which is much larger than the spontaneous emission rate into any given mode $\Gamma_{\downarrow}^{(m)}$ because the effective Purcell enhancement of spontaneous emission into any given cavity mode is weak. As the mode frequencies are below the molecular resonance we can approximate the critical molecular excitation as $f_{m}^{(c)}\gtrsim\Gamma_{\uparrow}^{(m)}/\Gamma_{\downarrow}^{(m)}$. Since, in addition, the molecules are not saturated, the molecular excitation can be approximated as $f_{m}\approx\Gamma_{\uparrow}/\Gamma_{\downarrow}$. Near threshold the molecular excitation must be close to the critical value, i.e. $f_{m}\approx f_{m}^{(c)}$, so that the inequality $\Gamma_{\uparrow}\Gamma_{\downarrow}^{(m)}\gtrsim\Gamma_{\downarrow}\Gamma_{\uparrow}^{(m)}$ follows. From the established condition $\Gamma_{\downarrow}\gg\Gamma_{\downarrow}^{(m)}$, one can thus conclude, that external pumping is far stronger than re-absorption of light from cavity modes, i.e. $\Gamma_{\uparrow}\gg\Gamma_{\uparrow}^{(m)}$. Below condensation threshold, where no mode population $n_{m}$ is macroscopically large, Eq. (3) can therefore be approximated as

[TABLE]

Within this approximation, Eq. (5) reduces to

[TABLE]

Equating this with $f_{m}^{(c)}$ (Eq. (7)) yields the threshold pump rate

[TABLE]

Given the assumption of sufficiently small mode populations $n_{m}$, Eq. (10) applies to the first condensation threshold only, i.e. it predicts the condensation of the mode $m$ with the lowest threshold pump rate. If this is the lowest mode $[0,0]$, then Eq. (10) describes the threshold of Bose-Einstein condensation; otherwise it predicts the onset of lasing.

Condensation into higher modes is more favorable than condensation into the ground mode, because the higher modes are closer to molecular resonance and therefore couple more strongly to the molecules. Thermalization of photons is therefore essential for the BEC phase as this redistributes photons from higher modes into the ground mode. If the system is condensed into a laser instead of a BEC, it is likely to remain a laser as pump rate increases, since the thermalization is unable to overcome the clamping of the molecules. This permits the estimation of the boundary between the laser and BEC phases for low pump powers, by considering the mode with the lowest threshold pump rate.

A.2 Boundary between single-mode and multi-mode condensation

The boundary between single-mode and multi-mode condensation can be derived very analogously to the discussion in Sec. A.1. Since, however, the population of the ground state mode is macroscopic, Eq. (3) can no longer be approximated by Eq. (8). Instead, it should be replaced by

[TABLE]

such that, Eq. (5) reduces to

[TABLE]

with the short hand notations $\Gamma_{\updownarrow}=\Gamma_{\downarrow}+\Gamma_{\uparrow}$ and $\Gamma_{\updownarrow}^{[0,0]}=\Gamma_{\uparrow}^{[0,0]}+\Gamma_{\downarrow}^{[0,0]}$.

Determining $f_{m}^{(s)}$ requires knowledge of the ground-state population $n_{[0,0]}$. The actual value of $f_{[0,0]}^{(s)}$ is close to its critical value, but does not exactly coincide with it, since $n_{[0,0]}$ is macroscopically large, but strictly speaking not diverging. Since Eq. (6) diverges for $f_{[0,0]}^{(s)}\to f_{[0,0]}^{(c)}$, it is not a good starting point for an estimate of $n_{[0,0]}$ for $f_{[0,0]}^{(s)}\lesssim f_{[0,0]}^{(c)}$. Eq. (12), on the other hand, can be solved for a finite value of $n_{[0,0]}$ by equating $f_{[0,0]}^{(s)}$ with $f_{[0,0]}^{(c)}$ as given in Eq. (7). The population of the ground mode below the second condensation threshold can thus be approximated by the relation

[TABLE]

with $f_{[0,0]}^{(s)}$ defined in Eq. (12). The corresponding value of $n_{[0,0]}$ then allows the calculation of the second condensation threshold by equating $f_{m}^{(s)}$ (for $m\neq[0,0]$) with $f^{(c)}_{m}$. The mode that achieves this threshold with the the lowest pump rate is then the mode that condenses in addition to the ground mode.

Evaluating these thresholds requires the explicit forms

[TABLE]

of the mode functions for a two-dimensional harmonic oscillator. The integrals can be evaluated analytically and read

[TABLE]

with

[TABLE]

and

[TABLE]

The $I_{k}$ can be expressed explicitly

[TABLE]

in terms of the polylogarithm functions

[TABLE]

The condition $f_{m}^{(s)}=f_{m}^{(c)}$ then determines the pump rate $\Gamma_{\uparrow}$ at which multi-mode condensation occurs. This condition can not be solved for $\Gamma_{\uparrow}$ analytically, but we found it to have a unique numerical solution.

With the principles of these analytic estimates, one can also explain which mode condenses at the lowest rate. The thermalization process redistributes photons from higher modes into lower modes, thus making it difficult for the higher mode to condense. At low thermalization rate mode $[0,1]$ condenses at the lowest pump rate to form phase ‘A’ in Fig. 2. As thermalization increases, molecules overlapping mode $[0,1]$ become clamped by the ground state, and any increase in the population of mode $[0,1]$ is transferred to the ground state. The thermalization between mode $[0,2]$ and $[0,0]$ is less than between $[0,1]$ and $[0,0]$, therefore increasing the thermalization coefficient has a smaller effect on $[0,2]$ than $[0,1]$. This allows mode $[0,2]$ to condense at a lower pump rate than $[0,1]$ for large thermalizations. This explains why the three-mode phase changes from ‘A’ to ‘B’ to ‘E’ as thermalization increases, and why a larger pump rate is required.

A.3 Decondensation

Last, but not least, let us discuss the decondensation of mode $[0,1]$ while modes $[0,0]$ and $[0,2]$ are condensed, i.e. the transition from ‘C’ to ‘B’ in Fig. 2.

Above the decondensation threshold (‘B’), modes $[0,0]$, $[0,2]$, and $[2,0]$ are condensed. We therefore replace Eq. (11) with

[TABLE]

where $S_{0,2}=\{\,[0,0],\,[0,2],\,[2,0]\,\}$. Eq. (5) therefore reduces to

[TABLE]

with the summations performed over $S_{0,2}$. As modes $[0,0]$, $[0,2]$ and $[2,0]$ are condensed their populations are determined from the condition $f^{(s)}_{m}=f^{(c)}_{m}$. This permits evaluation of $f_{[0,1]}^{(s)}$ following Eq. (15). Unlike the case of multi-mode condensation discussed above in Sec. A.2, the integrals can no longer be evaluated analytically, but Eq. (15) can readily be integrated numerically. In contrast to the previous cases, condensation is not found for increasing pump rate. It is rather obtained for decreasing pump rate, which implies decondensation with increasing pump rate.

Appendix B Absorption and Emission Spectra

All parameter values used for simulations match current photon BEC experiments 2017arXiv170706789D . The only difference between simulations and typical experiment values is in the mode spacing. We use a mode spacing of $30\,\mathrm{Trad/s}$ (i.e. $30/(2\pi)10^{12}\mbox{Hz}$) as this enables us to consider a system with few occupied modes, keeping it computationally tractable. Spacings in current experiments are a factor of $3$ smaller, but there is no fundamental difficulty in experiments with larger spacings.

The absorption (emission) of light from (into) a cavity mode with detuning $\delta$ by molecules is derived from experimental data nyman_robert_andrew_2017_569817 , shown in Fig. 3. The blue crosses indicate experimental data with a peak value of $1.2/s$. In our simulations we use $1.2\times 10^{-9}\kappa$ as peak value, and we use the fitted dependence satisfying the Kennard-Stepanov/McCumber relation, depicted in black.

The cavity modes have detunings of $\delta=\delta_{[0,0]}+m_{x}\omega_{x}+m_{y}\omega_{y}$ with $\omega_{x}=30\,\mathrm{Trad/s}$ and $\omega_{y}=30.3\,\mathrm{Trad/s}$, and $m_{x}$ and $m_{y}$ are the integers labeling each mode; $\delta_{[0,0]}$ is the detuning of the lowest mode from the molecular splitting and ranges from $-350\,\mathrm{Trad/s}$ to $-180\,\mathrm{Trad/s}$.

Appendix C Threshold

The threshold condition used in the simulations is an increase in the population of a mode greater than 3 orders of magnitude over an increase in pump rate of 10%. Here we justify our claim that the precise values used do not significantly change the phase diagram by presenting the large difference in mode population between the condensed and un-condensed modes.

Fig. 4 depicts a histogram for the occurrence of given mode populations. The data depicted in green corresponds to all modes considered un-condensed in the phase diagram depicted in Fig. 2, and data depicted in blue corresponds to modes categorized as condensed. As one can see, each condensed mode has a population larger than $10^{5.8}$ and each un-condensed mode has a population smaller than $10^{5.8}$, with the exception of 8 mis-categorized points with a population of about $10^{5}$. Inspection of these points reveals that they are part way through a condensation or de-condensation, and therefore constitute the shifting of a phase boundary by only a single sample. Changing the threshold condition to a change in population of $10^{3.4}$ instead of $10^{3}$ would result in perfect agreement between the threshold condition and the bi-modal distribution depicted in Fig. 4.

Appendix D Dependence on pump spot size

A natural parameter which is easily varied in an experiment is the width of the pump profile. In the main paper, we consider an essentially homogeneous pump profile. Narrow pump profiles have already been used to experimentally achieve multimode condensation, and they do so at lower total pump rates Marelic2016 since a larger fraction of the pump power is concentrated on the central, condensed modes. We therefore also investigated the dependence of the photon gas on the pump width, with a Gaussian pump profile centered around the center of the cavity.

If the width of the pump profile is narrower than the spatial extent of the lowest cavity mode, then pumping is mostly restricted to this mode. The majority of the pumped dye molecules are therefore clamped by this mode, which prevents multimode condensation. One therefore obtains single-mode BEC even for very strong pumping.

With increasing pump width excited modes also experience pumping, and sufficiently strong pumping results in multimode condensation. As the pump spot size is increased further, the pump power is distributed over more modes, and some modes lose their macroscopic occupation. There is thus an optimal pump spot size at which multi- mode condensation can be achieved with minimal pump rate. We found this size to be about $1.4$ times larger than the mean spatial extent of the ground-state wavefunction.