Photon-photon correlation of condensed light in a microcavity

The study of temporal coherence in a Bose-Einstein condensate of photons can be challenging, especially in the presence of correlations between the photonic modes. In this work, we use a microscopic, multimode model of photonic condensation inside a dye-filled microcavity and the quantum regression theorem, to derive an analytical expression for the equation of motion of the photon-photon correlation function. This allows us to derive the coherence time of the photonic modes and identify a nonmonotonic dependence of the temporal coherence of the condensed light with the cutoff frequency of the microcavity.


I. INTRODUCTION
A defining property of a Bose-Einstein condensate (BEC) is the macroscopic coherence exhibited by the particles in the lowest energy mode -a feature that allows the condensate to behave like a massive quantum wave [1,2].While features such as thermal equilibrium and large population of the ground state are the tell-tale signatures of a BEC, onset of quantum coherence can be a defining characteristic of condensates that do not thermalize completely and essentially operate out of equilibrium, such as exciton-polaritons [3,4] and photons [5][6][7][8].As such, theoretical and experimental investigation of coherence in both equilibrium and non-equilibrium condensates plays an important role in characterizing the properties of these macroscopic states.Over the years, coherence of BEC has been observed using interference experiments with ultracold atoms [9,10].Moreover, coherence in condensates of excitons-polaritons [11] and organic polaritons [12] have also been reported, while spontaneous phase selection [13] and spatio-temporal coherence [14,15] has been observed in photonic condensates.
The photonic Bose-Einstein condensates formed inside a dye-filled microcavity is driven-dissipative in nature, sustained by a detailed balance between the rate at which the dye molecules are incoherently driven and losses of molecular and photonic excitation in the system.The thermalization of the photon gas inside the microcavity is due to the vibrational states of the dye molecule, which allow thermal equilibration of photons via energydependent absorption and emission processes.From a physical perspective, one of the central differences between a photonic BEC and a laser is the rate of thermalization it is operating under.While a photonic BEC works in a near-equilibrium regime, close to thermal equilibrium with the dye molecules, a laser operates at a much lower thermalization rate and is firmly in the nonequilibrium regime.As such, macroscopic occupation of photons occurs in the lowest energy mode in a photon condensate, irrespective of the gain, whereas for a laser the large population corresponds to the mode with the highest gain.In between these regimes, lie exotic phases of light that exhibit strong multimode properties [16], non-stationary kinetics [17,18] and possible vortex-like features [19].
A key characteristic of Bose-Einstein condensates that operate in the quasiequilibrium regime is the spatiotemporal coherence.While condensates of excitonpolaritons exhibit significant inter-particle interactions, which lead to phase coherence [20] and superfluid behavior in the system [21], photons in a BEC do not interact strongly with each other [22,23].Instead coherence is created from indirect interaction of photons with the dye molecules, mediated via stimulated emissions [24].This leads to coherence properties such as symmetry-broken phase coherence [13], grand-canonical photon statistics [25] and transition from short-range to long-range spatial order across the condensation threshold [14,15] that have been experimentally observed.While phase coherence is generally well understood in these systems, studies on temporal coherence which is consistent with experimental results have been limited [14,16].
The properties of a photon condensation is well described by a microscopic model [26], derived from the light-matter interaction between a multimode cavity and an ensemble of emitters that represent the electronic and vibrational states of the dye molecules.The model can perfectly capture the thermalization [27] and quasiequilibrium properties [28] of the photon gas in the cavity, and also highlight features such as decondensation [29] and noncritical slowing down of dynamics [30].First order correlation function and the photon linewidth [27] based on a single-mode model predict a perfect temporal coherence in the condensed mode, consistent with the Schawlow-Townes limit at large photon numbers [31].While the theoretical results are in agreement with experimental observations in the vicinity of the BEC threshold [14,16], it exhibits significant deviation at higher photon numbers.
In this work, we derive the equation of motion of the photon-photon or the first order correlation function of the photon gas inside a dye-filled microcavity, especially when interactions between the different cavity modes are explicitly taken into account.Such a model already provides great insight into the role of spatial coherence in the kinetics of the condensed light in the cavity [19,28].Our focus here is on the temporal coherence of the condensed light, and by using the quantum regression theorem, we derive an analytical expression for the equation of motion of the photon-photon correlation.This allows us to study the temporal coherence or the coherence time of light inside the cavity for different system parameters.
The paper is arranged as follows.After the Introduction in Sec.I, we study the multimode model in Sec.II, and derive the rate equations for the photon correlations and molecular excitations.In Sec.III, we represent the photon-photon correlations in terms of the coefficients of the density matrix of the system.In Sec.IV, the time derivative of these coefficients are obtained to ultimately derive the equation of motion of the photon-photon correlation V.In Sec.VI, the time evolution of the correlation is numerically studied and the coherence time of the condensed light for different cutoff frequencies and pumping rate are presented.The results are discussed in Sec.VII.

II. THEORETICAL MODEL
The interaction of photons with the dye molecules inside a microcavity can be studied using a microscopic model [26].In its most general form, the theoretical model is valid for multiple cavity modes and also takes into account finite intermode correlations [28].The dynamics of the quantum system is governed by the following Master equation, where Ĥ0 = p δ p â † p âp is the bare energy of the cavity photons, and L[x] = 2xρx † − {x † x, ρ}, with {•} being the anticommutator.Here, âp (â † p ) is the annihilation (creation) operator of photons in the p th cavity mode, and κ is the rate at which it is lost from the cavity.The Pauli operators σ ± i denote the electronic states of the dye molecule at location r i in the cavity plane, which is pumped at a rate Γ i ↑ , but decays with a uniform rate Γ ↓ to non-cavity modes.The rate of absorption and emission of photons in the p th mode is given by A p and E p , respectively.The expression Ψ i p,q = ψ p (r i )ψ q (r i ) is given in terms of the mode functions ψ p (r i ) of cavity with quantum number p.The absorption and emission rates for the dye molecules, A p and E p , can either be calculated [27] or estimated from experimental data of the absorption crosssection [32].Moreover, these rates are known to follow the Kennard-Stepanov relation [33][34][35], A p = E p e −βδp , where δ p = ω ZPL − ω p , with ω p being the frequency of mode p and ω ZPL is the zero-phonon line.
The dynamics as well as the steady behavior of the cavity modes and the molecular excitation can be studied in terms of the equations of motion of the relevant observables such as the mode population or photon correlation from the above master equation.In most cases, it is easier to work with rate equations compared to directly solving the density matrix, especially when the system contains large number of modes and an inhomogeneous distribution of molecular excitations.
On the other hand, equations of motion for observables such as the mode population ⟨â † p âp ⟩ and the molecular excitation ⟨σ + i σ − i ⟩, or two-mode correlation function ⟨â † p âq ⟩ can be computed with relatively little computational resources when certain approximations are taken into consideration.First, it is helpful to use the fact that coupling between different dye-molecules have been neglected on account of the incessant collision between the molecules of the dye and the solvent, which quickly decoheres any intermolecule coherence.Secondly, one can invoke the semiclassical approximation such that correlations between molecules and photons can be factorized i.e., ⟨σ + i (t)â q (t)⟩ ≈ ⟨σ + i (t)⟩⟨â q (t)⟩.This is typically a good approximation when the number of emitters in a cavity is large, which is true in the case of dye-filled microcavity, where the number of dye molecules inside the microcavity is very large.
The photon population and correlations can be written as a matrix n, with elements n pq (t) = ⟨â † p (t)â q (t)⟩, and the molecular excitation fraction as a vector f with elements f i = j ⟨σ + j σ − j ⟩δ(r i − r j ).The semiclassical equations of motion for n and f are then given by, where f + and f − have elements A, E and Ω are diagonal matrices with elements A p , E p and ω p , respectively.Γ ↑ is simply a vector with terms Γ i ↑ .The matrices, Ẽ and Ã are diagonal with elements Tr[Φ i E(n+I)] and Tr[Φ i nA], respectively, where Φ i has the same dimension as n and has elements Ψ i pq .Importantly, Eqs. ( 2) and ( 3) allows for the computation of the photon correlation function ⟨â † p (t)â q (t)⟩ at time t, as well as at the steady state of the system, for a given set of system parameters.

III. PHOTON-PHOTON CORRELATION FUNCTION
The primary focus of the work is to calculate the firstorder correlation function, which would allow for the estimation of the spectral linewidth, as well as the temporal coherence of the multimode photon gas inside the dyefilled microcavity.The key quantity of interest is the photon-photon correlation, where p and q denotes the cavity modes.For most experiments involving the investigation of spatio-temporal coherence, the initial state of the system at t 1 is the steady state ρ ss .As such only the time difference, t = t 2 −t 1 is relevant and the two-time correlation reduces to c pq (t) = ⟨â † p (t)â q (0)⟩, where one can set the initial time t 1 = 0, and t 2 = t.Therefore, the two-time correlation can now be calculated using the term Tr[â † p (t)â q (0)ρ ss ].The correlation function can be calculated using the quantum regression theorem [37,38], which allows for setting up an equation of motion for the two-time correlation.The function ⟨â † p (t)â q (0)⟩ can considered as the time-evolution of the term ⟨â † p ⟩, governed by the master equation given in Eq. ( 1), however, with the initial state given by ρ ′ (0) = âq ρ ss .In other words, using the regression theorem, we have the relation, Similar approaches have been used to derive the first order [27] and second order correlation function [18].But these are mostly restricted to just single-mode systems, where there are no interactions arising from intermode correlations and the overall dynamics is fairly simple.
To apply the regression theorem, it is helpful to work with the density matrix formalism, which can be written in an expanded form using an appropriate orthonormal basis.For M photon modes and N molecules inside the cavity, one such basis is {|n, s⟩}, where |n⟩ = |n 0 n 1 . . .n p . . .n M ⟩ with n p being the population of the p th mode.Similarly, |s⟩ = |s 1 s 2 . . .s i . . .s N ⟩ is the molecular excitation state, where s i can take values 0 or 1, depending on whether the i th molecule is excited or not.Thus, the density matrix ρ can then be expressed in this basis as where |n, s⟩⟨n ′ , s ′ | = |n⟩⟨n ′ | ⊗ |s⟩⟨s|, based on the fact that there are no coherence between the molecules and the semiclassical approximation is valid.Now, the diagonal elements of ρ are those that satisfy n = n ′ and the off-diagonal terms are for n ̸ = n ′ .In general, n, n ′ ≥ 0, which means n p ≥ 0 ∀ p and implies that mode population is non-negative.As such, the following rearrangement of terms can be made, So, ρ can now be rearranged in terms of |n − k⟩⟨n|, where and R k n,s = ⟨n − k, s| ρ |n, s⟩.Note that k can be nonpositive and the above expression can be thought of as representing the density matrix by summing the k th diagonal terms.
The two-time correlation can now be computed by taking the trace Tr[â † p ρ ′ ], where where we expand ρ ′ similarly as in Eq. ( 8) with new set of coefficients {P k n,s } which are related to the steady state coefficients by where k q is defined as a vector with the q th term as unity and zero elsewhere.As such, we obtain the following expression, This gives us a time-evolution of the photon-photon correlation in terms of the coefficients,

IV. CALCULATION OF COEFFICIENTS
The estimation of the photon-photon correlation now depends on finding the solution to Eq. (11), which is obtained by finding the terms Ṗ k p n,s .In particular, the above time derivative needs to be expressed as a function of the coefficients defined in Eqs. ( 8) and ( 9), which gives the expression for the operator ρ ′ in terms of {P k p n,s }.Hence, to obtain the time-derivative of the two-time correlation, the natural step is to calculate the derivative of the operator ρ ′ using the master equation in Eq. ( 1).In particular, the different terms in the master equation needs to be investigated for their contribution to the time-derivative, starting from the Hamiltonian Ĥ0 to the Lindblad terms L[σ ± i ] and L[â p ].For example, contribution from Ĥ0 in finding the relevant terms in ρ will simply come from where Ĥ0 = m â † m âm .Here are the contributions from the different terms in the master equation: • where C 0 (m) = (n m + 1)(n m + 1 − δ p,m ).Note that κ describes the loss rate for all cavity modes, and a state n can lose a photon in mode m and thus change into n − k m , where k m (as defined earlier) is a vector with the m th element as 1 and 0 everywhere else.Conversely a photon population n can come from state n + k m by losing a photon in mode m.As such, these states have probability P where Γ↑ = i,∀si=0 Γ i ↑ and s i is a vector with the i th element as 1 and 0 everywhere else.As Γ i ↑ is the pumping rate it excites the state s to s + s i by exciting the i th molecule.Similarly, a state s may be pumped from a state s − s i with probability being decay of molecule with rate Γ ↑ , simply produces a the reverse state transformation.
A more detailed description of the above transformations is shown in Appendix A. where Note that the process of absorption and emission, as governed by A m and E m ′ in Eq. ( 1), gives rise to intermode correlations.In the case of m = m ′ the mode is coupled to itself.Mathematically, this is a state with probability P k p n,s , which upon absorption of a photon in mode m by the i th molecule at location r i , will transform to a state with probability P k p n+k m ,s−s i .Now, including m ̸ = m ′ , introduces inter-mode correlations, which reveals itself in the coupling of different k i terms.For instance, k p in the derivative on the left hand side is connected to the k p −k m ′ +k m on the right for all m and m ′ .Similar calculations can also be done for transformations arising from the emission term E m ′ .The contribution from Hermitian conjugates in Eq.1 are calculated in similar manner.Now, the full expression of Ṗ k p n,s is simply the sum of the different contributions, which gives us a general relation of the time derivative to the set {P k p n,s }.

V. EQUATION OF MOTION
In Eq. (11), we represent the time-derivative of the two-time correlation function, ċpq (t) = d dt ⟨â † p (t)â q (0)⟩ in terms of the rate of change of the probabilities Ṗ k p n,s .Moreover, in Eqs. ( 13)-( 16), we derived the timederivatives in terms of the probabilities {P k p n,s }.As such, the equation of motion of the two-time correlation function ⟨â † p (t)â q (0)⟩ can now be derived independently of these probabilities, which in practical conditions can be very difficult to estimate.
From Sec.IV, the contributions from the Ĥ0 and L[â p ] terms in the master equation lead to and therefore introduce an oscillatory and a decay term in the equation of motion.The terms related to pumping and decay of molecules in the system, given by L[σ ± i ], do not contribute to the photon correlation.However, the absorption and emission terms, given by A m and E m ′ , make a significant contribution to the equation of motion.For A m this is given by The contribution from E m ′ is given simply by replacing absorption term with emission, the index m by m ′ , and the summation over all s i = 1.An extended derivation of Eq. ( 18) is shown in Appendix B.
Importantly, the contributions from A m and E m ′ still contains coefficients P k m n,s .A useful condition to use at this point is the semiclassical approximation discussed in Sec.II, which factorizes the correlation between photons and molecules.This leads to P k p n,s ≈ P k p n P s , where one can interpret P s to be probability for molecules to have an excitation profile described by s.As such Eq.( 10) can be written as, where the sum over probabilities of all possible molecular excitation profile is unity, i.e., s P s = 1.Now, to simplify the expression in Eq. ( 18), we use the semiclassical approximation in Eq. ( 19), such that ċpq The coefficient s (P s i,∀si=0 Ψ i p,m ), here is not straightforward.It is a sum of the mode function Ψ i p,m at locations r i , where the molecule is not excited as described by the vector |s⟩.The above expression can be rewritten as shown in the Appendix C, as i Ψ i p,m (1 − f i ), where f i is the probability of finding an excited molecule at location r i or the excitation fraction.A similar expression can also be obtained from the emission E m ′ term.
Therefore, bringing all the terms together, the full equation of motion of for the photon-photon correlation is given by where The photon-photon correlation or the first order correlation function allows for the estimation of the emission spectrum of the cavity and the temporal coherence of the different modes.The time evolution of the correlation function c pq (t) = ⟨â † p (t)â q (0)⟩ is governed by the equation of motion in Eq. ( 21), which can be numerically solved for a particular set of system parameters.The initial state of the system, at time t = 0, is the steady state correlation given by ⟨â † p (0)â q (0)⟩ = ⟨â † p âq ⟩ ss , which is obtained by finding the the steady state solutions of Eqs. ( 2) and (3) discussed in Sec.II.

VI. TEMPORAL COHERENCE
The temporal coherence of the different photonic modes can be studied from the spectral function of the emitted light from the cavity [36], given by the relation, which can be estimated by numerically solving for ⟨â † p (t)â p (0)⟩ using Eq. ( 21). Figure 1 shows the emission spectrum of the condensed ground state mode p = 0, in comparison with the first and second excited modes, p = 1 and p = 2 respectively.The plots show that the ground state has a much narrower linewidth compared to the higher modes, which implies high temporal coherence in the condensed ground state mode.A notable point is the spectrum of the excited modes.The first  0)⟩ = ⟨np(0)⟩ at t = 0), for the ground state mode (p = 0) and the first and second modes (p = 1, 2).The pump is considered to be a Gaussian focused at the center, with a width equal to thrice the oscillator length l0.The cavity cutoff frequency ω0 ≈ 520 THz, with spacing between two adjacent cavity modes, ∆ω = 1.7 THz.The rate of loss of cavity photons is equal to κ ≈ 0.2 THz and the absorption (Am) and emission (Em) rates are calculated from experimental data [32].The loss of excitations outside the cavity modes is Γ ↓ ≈ 3 × 10 −5 THz and pumping rate Γ ↑ = 0.2 × Γ ↓ .
2. Time evolution of the first order correlation function.
The figure shows the photon-photon correlation of the lowest energy state (p = 0) as a function of time t, normalized with the initial steady state population ⟨â † p (0)âp(0)⟩ = ⟨np(0)⟩ at t = 0.The plots show the temporal behavior of the function for three different pump rates Γ ↑ (in units of Γ ↓ ).The bold circles mark the coherence time τ0 for the different pump rates.All other system parameters are the same as in Fig. 1.
excited mode has lower emission than the second excited mode, which implies that the second mode is more populated.This is due to the choice of a Gaussian pump spot focused at the center of the cavity, which tends to excite the even modes, with mode-functions that peak at the center.Moreover, the even modes also couple more strongly to the ground state [39].Coherence time τ 0 (units of ps) 3. Variation of temporal coherence with cavity cutoff frequency.The plots show the coherence time of the lowest energy or ground state mode for different cutoff frequencies of the cavity, and for three different pumping rates Γ ↑ (in units of Γ ↓ ).All other system parameters are the same as in Fig. 1.
Figure 2, shows the time evolution of the photonphoton correlation of the ground state (p = 0) for different pump powers.To facilitate the comparison, the correlation function is normalized with the initial steady state population at t = 0. photon-photon correlation decreases with time, but the rate of decrease reduces as the pump power is increased, which confirms the increased temporal coherence of the photons at higher pump powers.Importantly, from the close to Lorentzian lineshape in Fig. 1, it can be inferred that the decrease of the correlation function is exponential.Thus the coherence time τ p for cavity mode p can be defined from the photon-photon correlation as The temporal coherence τ 0 for the ground state for different pump powers are shown in Fig. 2 (bold circles), which marks the time at which the correlation function c p (t) drops to 1/e times its initial value c p (0). Figure 3, shows the variation of the coherence time τ 0 of the lowest energy state with the cutoff frequency ω 0 of the cavity.The cutoff frequency of the cavity is a critical system parameter that not only defines the energy of the lowest cavity mode, but also photon-energy dependent rates of absorption (A p ) and emission (E p ) of photons by the dye-molecules inside the cavity.The figure shows that the coherence time is dependent on the cavity cutoff frequency, but does not vary monotonically with change in the cutoff frequency.This is either due to the nonmonotonic variation of thermalization inside the cavity across the chosen range of cutoff frequencies or strong mode competition for molecular excitations.The figure also shows how the coherence time of the ground state depends on the pump power, which controls the total number of photons in the cavity and therefore drives the photon condensation transition.

VII. CONCLUSION
In this work, we derive an equation of motion for the first order correlation function or the photon-photon correlation for the photon gas inside a dye-filled microcavity.The nonequilibrium model takes into account a multimode photonic cavity, where finite inter-mode coherences are not completely ignored, which makes the calculations significantly more complex, but allows us to compute photon-photon correlations between different modes.Importantly, these relations allow the theoretical and computational investigation of temporal coherence of photon condensates that are consistent with actual experimental findings for a far wider set of parameters [39].
The work opens the door to study more complex behavior of photon correlations and temporal coherence, especially in regimes where strong mode competition for excitation exists.A particular phenomenon of interest is to study how temporal coherence behaves in the regime of multimode condensate, where molecular excitations are clamped by a condensing mode and different modes compete to unclamp the excitation, leading to the phenomenon of decondesation [29].Moreover, other directions include a more comprehensive study of spatiotemporal correlations in photon condensation under the influence of non-stationary pump, where phenomena such as vortex-like structure formation [19], and partially coherent light can be engineered.pumping and decay of molecules as governed by the rates Γ i ↑ and Γ ↓ , respectively.These are given by n,s The first term on the right hand side (RHS) of Eq. (A1) contains the contribution from the {σ − i σ + i , ρ} term of the operator L[σ + i ] which acts on |s⟩⟨s|, if and only if s i = 0 i.e., σ − i σ + i |s⟩⟨s| = |s⟩⟨s|σ − i σ + i = δ si,0 |s⟩⟨s|.Comparing the basis states |n − k p ⟩⟨n| and |s⟩⟨s| on both sides, the coefficient from the first term is −P k p n,s Γ↑ , where the term Γ↑ = i,∀si=0 Γ i ↑ , which is the total pump rate at all unexcited sites in the state |s⟩.The second term on the RHS arises due to the contribution from the σ + i ρσ − i term of L[σ + i ].Note that this term increases the excitation at locations where the molecules are unexcited or σ + i |s⟩⟨s|σ − i = δ si,0 |s + s i ⟩⟨s + s i |.This is illustrated in Fig. 4. As such we obtain the term i,∀si=0 Γ i ↑ |s + s i ⟩⟨s + s i |, which for any fixed s is a summation over all states with one more excitation than s.Taking all s into account and noting that there will be no |s⟩⟨s|, where s i = 0 ∀i, the total contribution of second term maybe rewritten by a simple change of indices n,s The change of indices above is better illustrated in Fig. 5. Now by comparing terms in LHS and RHS with the same basis |n − k p ⟩⟨n| ⊗ |s⟩⟨s|, we obtain The contribution from decay to non-cavity modes governed by the Γ ↓ term can be obtained in a similar manner.
Next, there are the contributions due to the absorption and emission processes, which are given by operators of the respectively.For the A m term, the relevant relation between the coefficients is n,s Calculating the first term in RHS gives us the contribution from the term âm where n − k m = n ′ and the indices in the molecular states have been rearranged.
Similarly, the second term in RHS is the contribution due to â † Combining the relations from Eqs. (A5) and (A6), align the corresponding basis, the following expression for the coefficient is obtained The contribution from the Hermitian conjugate and molecular emission, given by the rate E m ′ in the master equation, can be calculated along similar lines.
as the m ̸ = l terms cancels out in the calculation.Next, the contributions from the terms L[σ ± i ] are considered: For each s ′ on the right there exists only one term from left with s = s ′ + s i (∀s i = 0), which then gives c s ′ = i,∀si=0 Γ i ↑ .Including the expression for c s ′ above, the following relation is obtained: A similar derivation can be done for contribution from the term L[σ − i ].The contribution for the absorption and emission, given by the rates A m and E m ′ , respectively, are considered next.First, the term A m for the case m = m ′ : where |ψ i m | 2 is squared wavefunction of cavity mode m.Now, the last two lines in Eq. (B6) are all m ̸ = p contributions, and by taking n ′ = n + k m , the last two lines become  For m ̸ = m ′ , there can be two cases, either m = p or m ̸ = p for any m ′ .For the latter case, d dt ⟨â † p (t)â q (0)⟩ = 0.For m = p, the expression is Again, similar calculations exist for terms arising from emission, E m .
Applying the approximation in Eq. (C1) to the equation of motion: d dt ⟨â † p (t)â q (0)⟩ = (iδ p − κ 2 )⟨â † p (t)â q (0)⟩ − The term s (P s i,∀si=1 Ψ i m,p ) can be transformed into i Ψ i m,p ( s ′ ,∀s i =1 P s ′ ), where the last summation is over all s ′ with unity at position i.The sum s ′ P s ′ is simply the total probability of net excitation of the i th molecule (i.e. the molecules at position r i ).Using a similar argument, s (P s i,∀si=0 Ψ i p,m ) can be changed to i Ψ i p,m ( s,∀si=0 P s ), where s P s is the net probability of i th molecule being unexcited.Let f be a vector, where f i is the excitation fraction at position r i , which results in i Ψ i p,m ( s ′ P s ′ ) = i f i Ψ i p,m and i Ψ i p,m ( s P s ) = i (1 The equation of motion of the photon-photon correlation under the semiclassical approximation is then given by,

2 FIG. 1 .
FIG.1.The emission spectrum of the first few photonic modes.The figure shows the spectral function Sp(ω), normalized with the initial steady state population ⟨â † p (0)âp(0)⟩ = ⟨np(0)⟩ at t = 0), for the ground state mode (p = 0) and the first and second modes (p = 1, 2).The pump is considered to be a Gaussian focused at the center, with a width equal to thrice the oscillator length l0.The cavity cutoff frequency ω0 ≈ 520 THz, with spacing between two adjacent cavity modes, ∆ω = 1.7 THz.The rate of loss of cavity photons is equal to κ ≈ 0.2 THz and the absorption (Am) and emission (Em) rates are calculated from experimental data[32].The loss of excitations outside the cavity modes is Γ ↓ ≈ 3 × 10 −5 THz and pumping rate Γ ↑ = 0.2 × Γ ↓ .

FIG. 4 .
FIG. 4. Visualization of the term i,∀si =0 Γ i↑ as a square, where each block ri denotes a position, and black (white) implying the molecule at the position is excited (non-excited).For a specific s there are associated squares with one more excitation, give by a term Γ j ↑ which creates excitation at a new position rj, with new basis |s + s j ⟩⟨s + s j |.

FIG. 5 .
FIG. 5. Visualization of |s ′ ⟩⟨s ′ | as a square, with different contributing terms coming from a set of squares, with one less excitation at position rj compared to s ′ .Thus, each term contributes Γ j ↑ P k p n,s ′ −s j .