Phase transitions of light in a dye-filled microcavity: observations and simulations

Photon thermalisation and condensation in dye-filled microcavities is a growing area of scientific interest, at the intersection of photonics, quantum optics and statistical physics. We give here a short introduction to the topic, together with an explanation of some of our more important recent results. A key result across several projects is that we have a model based on a detailed physical description which has been used to accurately describe experimental observations. We present a new open-source package in Python called PyPBEC which implements this model. The aim is to enable the reader to readily simulate and explore the physics of photon condensates themselves, so this article also includes a working example code which can be downloaded from the GitHub repository.


Photon condensation in dye-filled microcavities
Light in resonators which are filled with a fluorescent medium has been studied in many guises, with the most notable variant being the laser. More recently attention has turned to what happens when the medium also interacts strongly with the light, especially when the resonator supports multiple modes.
When the coupling between light and matter is coherent, this is known as the strong-coupling regime of cavity quantum electrodynamics (CQED). The resulting hybrid light-matter particles are known as polaritons, and if pumped sufficiently strongly can cause lasing [1,2,3]. The particles may thermalise among themselves to form polariton condensates, which has been shown in planar dielectric microcavities filled with inorganic semiconductors [4] or organic dyes [5,6], and also plasmon-polariton metamaterial systems [7].
There is a simpler, incoherent form of strong coupling between light and matter, which is reached when the absorption of light by the gain medium is at least as fast as loss from the resonator. This is the regime in which photon thermalisation to room temperature and Bose-Einstein condensation (BEC) [8] can occur, and this regime is the topic of this work. The basic principles have been described elsewhere, but we will start by re-stating them in a more general form. We then review some of our group's recent experimental and theoretical results, with the perspective that photon condensation is an example of a wider class of dynamical phase transitions, from which one can learn a lot about phase transitions in general. Finally, we present a new suite of open-source software that we have published to simulate photon condensation called PyPBEC, including a result on the statistical properties of multi-mode microlasers derived using the software.

How light reaches thermal equilibrium
Optical resonators (and metamaterials) modify the density of states (DoS) available for light. We consider here a class of resonators which provide a large gap in the DoS, with several modes lying just at energies above the gap energy. Furthermore, the DoS increases with energy, such that there is exactly one lowest-energy state if the the system is finite in size. The simplest such resonators are open microcavities of length a few half-wavelengths of the light and radius of curvature at least a few tens of microns.
A fluorescent medium is introduced to interact with the light. The medium is chosen so that the absorption and emission spectra do not extend beyond the gap in the DoS for light. Thus the combination of resonator and medium provide a robust energy minimum for the light. They also overlap slightly, in that there is an energy range where both absorption and emission of light can occur rapidly. In that overlap region, we require that the ratio of absorption and emission spectra follows a Boltzmann factor: absorption dominates at higher energies, emission dominates at lower energies. This relation is valid for many fluorescent dyes (known as the Kennard-Stepanov relation), some doped laser glasses (known as the McCumber relation) and many semiconductors (a limit of the van Roosbroeck-Shockley relation).
If the resonator lifetime is longer than the absorption timescale, and most of the fluorescence light is emitted into the resonator, the average population of photons on the resonator modes depends on the ratio of absorption to emission. Since we chose the material for the Boltzmanntype relation between absorption and emission, the photon populations will follow a thermal factor. When the populations increase, stimulated emission means that the Bose-Einstein distribution is the valid distribution (i.e. with population controlled by the chemical potential). With sufficient population, the number in the lowest-energy state above the gap energy will diverge, while numbers of excited states saturate. Thus, Bose-Einstein condensation is observed.  (Right) Experimental data on the occupation of low-lying energy modes in a small photonic system (data re-purposed from Ref. [9]). Beyond a critical pump power, a significant fraction of the total photons condense to the ground state.
Several groups have achieved photon thermalisation and BEC using dye-filled open microcavities [8,10,11,12] as shown in Fig. 1. An open microcavity consists of a pair of high-reflectivity, dielectric mirrors spaced a few half-wavelengths of light apart, at least one of which is curved. Light propagating has, in the paraxial approximation, a density of states exactly as described above, specifically equivalent to a two-dimensional harmonic oscillator (with DoS increasing linearly with energy above the ground state). The cavity is filled with a fluorescent dye such as Rhodamine 6G which is optically pumped. With a ground state energy of 2.1-2.2 electron-Volts (equivalent to 560-590 nm) the absorption by the dye is faster than cavity loss, so thermalisation occurs. For sufficiently strong pumping, BEC of photons is observed.
More recently, some other physical implementations of photon thermalisation and condensation have been developed. Plasmonic lattices can be used to enhance light-matter couplings with band-structures and materials suitable for BEC, as demonstrated by Hakala et al [13]. Dopings in fibre lasers also have the appropriate relation between the absorption and emission spectra, and absorption is strong, so BEC is possible with appropriate intra-cavity filtering [14]. Semiconductors built into microcavities are commonplace, and some evidence has come to light that pre-existing planar microcavity devices may show thermalisation and BEC [15].
As this article is intended as a conference proceedings, rather than reviewing the complete literature, we will now specialise to our own recent results, explaining the meaning rather than the details.

Phase transitions in steady state
Bose-Einstein condensation is an example of a phase transition. In the thermodynamic limit of a great many particles at thermal equilibrium, there are discontinuities in several properties as a function of the thermodynamic parameters (temperature and chemical potential). Photons thermalising in dye-filled microcavities constitute a finite-sized, driven-dissipative system which can, at best, only approximate thermal equilibrium. Therefore, not only are the phase transitions modified from the simple BEC case, but perhaps even the concept of phase transition should be re-evaluated.

Breakdown of thermal equilibrium, including decondensation
The thermalisation rate is effectively the rate at which photons are absorbed by the dye. This rate can be compared to the cavity photon loss rate to judge how nearly the system reaches thermal equilibrium. We made a numerical and theoretical study of the non-equilibrium phase diagram of a pumped dye-filled microcavity [16]. Steady-state photon populations were categorised as functions of pump and absorption rates relative to cavity loss rate. For the simulation parameters used, the population of a mode was found to jump by at least a factor 100 in a range of about 10% variation in pump rate, which provided a working definition of a phase transition.
When absorption was at least 10 times faster than loss, the only phase transition detected was BEC. The opposite case (extremely weak absorption) lead to a series of phase transitions corresponding to single-or multi-mode lasing. Which modes transitioned depended sensitively on control conditions (pump spot size and position, thermalisation rate, pump rate).
The intermediate, partially-thermalised case (with absorption and cavity loss rates being roughly equal) is perhaps the most interesting. As pump rate is increased, the first mode to condense (jump up in population) is robustly the ground state, but then further modes may condense. For very large pump rates, competition between modes for the gain given by the molecular excited state population leads to an effect we called decondensation. In that case, the population of a mode decreased sharply with increasing pump rate.

Fuzzy phase diagrams for small systems
We have subsequently explored the same phase diagram experimentally [9]: see Fig. 2. Our apparatus allows us to measure the populations of individual energy levels, and we use those populations to assign phases to regions of experimental-parameter space (pump laser power and cavity length).
Unlike the simulations, the populations are small even at the threshold for condensation or lasing, with as few as 7 photons required for BEC [17]. The phase transition is no longer a clear discontinuity in an order parameter, but rather is evidenced by slightly super-linear variation in population with respect to pump power. The phase transition cannot be assigned using heuristics (informed guess-work) so we turn to machine learning.
The populations over the parameter space show clustering, and we apply a fuzzy c-means algorithm to identify those clusters (the colours in Fig. 2). Furthermore, we use a measure we call "membership entropy" to quantify the degree of affinity between population configurations and phases, i.e. how strongly we can say that a configuration sits within a phase. This measure is shown by the level of whiteness in the figure and represents the ambiguity in assigning configurations to phases. In this way, we construct a fuzzy phase diagram.

The dynamics of phase transitions
We now turn to dynamical results (away from steady state), after quenches (theory work) and pulses (experiments). Close to phase transitions, the time taken to reach the steadystate behaviour increases dramatically, a phenomenon known as "critical slowing down". In simulations we observe exactly that behaviour [18] for most of the BEC and multi-mode condensation transitions, except the decondensation. While the decondensation is associated with critical slowing down close to threshold, far above the threshold pump rate, there persists a slowing down by at least an order of magnitude with respect to the expected timescales of the system. This non-critical slowing is related to the strong gain competition which causes the decondensation. The dynamics in photon condensation normally happens in 5-500 ps. It is not feasible to quench a pump power between two values that fast, so we use 40-ps pulses to drive dynamics [9]. We then observe the dynamics using single-photon detectors (although operating in a regime of many photons) and time-to-digital conversion electronics which have timing resolution as short as 40 ps. First, we observe the pulse-averaged behaviour to observe the total populations; we detect threshold behaviour as in the steady-state case. Then, we measure the average population as a function of time after the pulse (see Fig. 3), and we detect the transient analogue of critical slowing down. Lastly, we use a non-stationary, two-time correlation measure g (2) (t 1 , t 2 ). From the anti-correlations between early and late photon events, we note that the slowing down is related to strong jitter. The experiments match theory rather well. We conclude that the BEC phase transition is caused by (deterministic) stimulated scattering/emission but can only occur after (stochastic) spontaneous scattering of photons into the condensing mode.

Generalised simulation of light in multi-mode resonators with fluorescent media
Photon condensation is a rich phenomenon which is driven by both experiments and simulations. As such, it is important to have an accurate theory of the system based on a fundamentally We have such a model, which successfully predicts and explains many of our experimental observations [9,17,19,20]. We have shared our software for these simulations as an open-source project called PyPBEC (Python P hoton B ose-E instein C ondensation, available at https://github.com/photonbec/PyPBEC). Despite the name, it can be applied to a large range of physics involving multi-mode optical resonators interacting with fluorescent media, such as microlasers or plasmonic lattices.

The model used and its limitations
The model originates from work by Kirton and Keeling, where they start from a multi-mode, multi-molecule CQED framework, including the effects of mechanical relaxation (phonons) associated with the molecules [21,22]. They extended it to include effects of spatial inhomogeneities in the excitation of molecules [23]. The model is defined by the equations: whereρ is the density operator of the photon-molecule system. Here,Ĥ 0 is the energy-conserving Hamiltonian and L[x]ρ =xρx † − 1 2 x †x ,ρ is the standard Lindblad operator. Moreover,â m andâ † m are the photon annihilation and creation operators for the cavity mode m, whereas σ ± l are the lowering (-) and raising (+) Pauli operators for the lth molecule. A m and E m are the absorption and emission rates for molecules into or from mode m, and Γ ↑ is the rate of incoherent pumping of the molecules, κ is the rate of photon loss and Γ ↓ is loss rate due to fluorescence in all non-cavity modes or other non-radiative processes.
In the form we use the equations most often (with a mean-field, semiclassical approximation assuming molecule and photon populations are uncorrelated), the above master equation can be used to derive a set of coupled equations of motion for the photon and molecular excitations: The molecules in the cavity have been divided into spatial bins j, with f j being the fraction of molecular excitations at j. The photon population in mode m is given by n m . Furthermore, g m,j is the coupling between the mode m with molecules in the spatial bin j. We note that in the above equations we assume that the mode-mode coherence is weak and can be ignored if one is only interested in studying the photon and molecular excitations.

The open-source code and its structure
The code for simulating equations of motion is structured in three components: • The photonic modes: the Cavity superclass defines all the energies and spatial representations of the resonator modes. It also defines the representation to be used if a spatially-resolved simulation is being used.
• The fluorescent medium: the OpticalMedium class facilitates calculation of values for absorption and emission spectra at requested energies. A subclass of OpticalMedium called Rhodamine6G for obvious reasons is implemented, based on open-source data at [24]. Emission and absorption rates can alternatively be input directly.
• The solver (which defines and solves the equations of motion): the Solver superclass defines a format for recording the equations of motion, initial conditions, and the computational methods for solving the equations. Initially, three subclasses are defined, SteadyState for solving the equations above for the steady-state photon and molecule populations and ODE for the dynamics, as well as a MonteCarlo class which solves stochastically for the populations of the master equation without making the mean-field approximation.
The packages are written in Python3. A typical program starts by importing the three components. The OpticalMedium is either selected from the available list, custom defined, or emission and absorption rates at selected wavelengths may be input directly to the Cavity object. Once the Cavity is defined, the cavity parameters such as κ and g l,m are set, including the rates of pumping Γ ↑ . The Solver object is instantiated, and assigned a Cavity object and initial conditions. Finally the Solver.solve method is invoked, and results can be analysed and visualized.

Notes on extensibility
The code is structured so that the medium, the photonic modes, the coupling between medium and photonic modes, the equations of motion and the solution method can be easily customised. As such, it is extremely general. The specialisation to microcavity problems and photon BEC in particular comes about through the examples given, and the solvers, cavities and media which have been implemented.
At present, we have only implemented Rhodamine 6G as a medium, but other fluorescent dyes media based on many two-level emitters which undergo rapid dephasing can easily be added by including experimental data or simple theoretical models. The cavities can be any optical structure which exhibits discrete modes.
The choice of solver will of course play a significant role in both the results which can be obtained and the computational cost. The MonteCarlo solver runs much slower than the SteadyState solver, but gives beyond-mean-field results which match well to data, see Ref. [19].
We plan to implement further solvers which take account of inter-mode phase coherences or use quantum-regression methods to efficiently calculate beyond-mean-field effects. Likewise, further materials and cavity types are also planned.

Example simulation: Stochastic dynamics of a two-mode laser
In this section, we give an example of how PyPBEC can simulate the stochastic dynamics of a pair of degenerate lasing modes inside a dye-filled resonator. We use the cavity lifetime 1/κ as the unit of time, and consider degenerate absorption and emission rates, E m = 50κ and A m = 10κ for m = 0, 1, with molecular pump and decay rates Γ ↑ , Γ ↓ = {100κ, κ}. The mean-field evolution of the two degenerate modes is obtained using the ODE solver to obtain the steady state of the system (Fig. 4, left panel). The stochastic dynamics starting from the steady state are obtained with the MonteCarlo solver. Fig. 4, right panel, shows how the photonic excitations switch between the two degenerate modes during one such stochastic evolution; populations are anticorrelated. The code to generate the steady-state and stochastic dynamics shown in Fig. 4, is given in the example folder of the PyPBEC package. The high-level code is shown here. First, the Cavity class object is defined (M and J are the numbers of cavity modes and spatial bins respectively). import numpy as np from PyPBEC.Cavity import Cavity cavity = Cavity(M=2, J=1, cavity_loss_rates= [1,1] An OpticalMedium object is not needed as we have explicitly defined absorption and emission rates. Then we solve the ordinary differential equation (ODE) to find mean-field evolution towards the steady state: It is then a simple matter to plot the results, using the .t and .photons attributes of the resulting Solver objects.