## Fluctuation dynamics of an open photon Bose-Einstein condensate

Fahri Emre Ozturk, 1, ∗ Tim Lappe, 2, † G¨ oran Hellmann, 1 Julian Schmitt, 1, ‡ Jan Klaers, 1, § Frank Vewinger, 1 Johann Kroha, 2 and Martin Weitz 1 1 Institut f¨r u Angewandte Physik, Universit¨t Bonn, Wegelerstr. a 8, 53115 Bonn, Germany 2 Physikalisches Institut and Bethe Center for Theoretical Physics, Universit¨t a Bonn, Nussallee 12, 53115 Bonn, Germany

(Dated: October 9, 2019)

Bosonic gases coupled to a particle reservoir have proven to support a regime of operation where Bose-Einstein condensation coexists with unusually large particle-number fluctuations. Experimentally, this situation has been realized with two-dimensional photon gases in a dye-filled optical microcavity. Here, we investigate theoretically and experimentally the open-system dynamics of a grand canonical Bose-Einstein condensate of photons. We identify a regime with temporal oscillations of the second-order coherence function g (2) ( τ ), even though the energy spectrum closely matches the predictions for an equilibrium Bose-Einstein distribution and the system is operated deeply in the regime of weak light-matter coupling. The observed temporal oscillations are attributed to the nonlinear, weakly driven-dissipative nature of the system which leads to time-reversal symmetry breaking.

## I. INTRODUCTION

Bose-Einstein condensates are the experimental basis of a variety of observed collective quantum phenomena [1, 2]. Generally, Bose-Einstein condensation is a phenomenon in thermodynamic equilibrium for Bose systems, usually described for a fixed total particle number, which leads to a macroscopic ground-state occupation. In ultracold atomic gases [3, 4] and exciton-polariton systems [5, 6], condensation has been observed following thermalization by interparticle collisions, a process that leaves the total particle number constant. Photons usually do not exhibit condensation, with the chemical potential vanishing at all temperatures in the wellknown example of a (three-dimensional) blackbody radiator. Two-dimensional photon gases under harmonic confinement can, however, reach condensation, as was experimentally demonstrated in dye-filled optical microcavities [7-9]. In this system, thermalization is reached by absorption and re-emission processes on the dye molecules, which leaves the average particle number constant but allows for fluctuations around this average value. In the limit of number fluctuations that become as large as the average particle number, such a situation is described well by the grand canonical statistical ensemble [10-13]. Lasers are long-known physical systems that also have macroscopic population of excited states, but which operate far from thermal equilibrium [14, 15]. The crossover from lasing to condensation, characterized by a varying

∗ Email: oeztuerk@iap.uni-bonn.de

† Email: lappet@th.physik.uni-bonn.de

‡ Present address: Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom § Present address: Complex Photonic Systems (COPS), MESA+ Institute for Nanotechnology, University of Twente, 7522 NB Enschede, The Netherlands

degree of thermalization, has been investigated by observing deviations from a thermalized distribution, both in polariton and photon gases [8, 16-20]. Particle-number conserving condensates as well as lasers are characterized by vanishing particle-number fluctuations in the thermodynamic limit, i.e., by a value of the second-order coherence function of g (2) ( τ ) = 1 at all delay times τ .

Experimentally, evidence for Bose-Einstein condensation in the grand canonical statistical regime with unusually large statistical number fluctuations has been observed with photons in the dye microcavity system [21]. The coupling of photons to the photo-excitable dye molecules implies that the dye does not only act as a heat bath, but also as a particle reservoir due to the possible interconversion of photons and dye electronic excitations. For a large relative size of the dye reservoir, this leads to strikingly enhanced statistical number fluctuations and a zero-delay second-order correlation g (2) (0) = 2, i.e., as in a thermal source. Notably, these fluctuations, which can be as large as the average value, occur deep in the condensed phase. On the other hand, with a smaller effective relative size of the dye reservoir, the dye microcavity photon condensate can also be operated in the (usual) canonical statistical regime, with much smaller number fluctuations and a zero-delay intensity correlation g (2) (0) = 1. Independent of the ensemble conditions, frequent collisions of solvent molecules with the dye on a timescale of 10 -14 s [22], cause the dye microcavity condensate to operate in the weakly coupled regime of matter and light [23, 24], i.e., the trapped particles are photons (not polaritons) and the system can be well described by a rate equation model, with, e.g., no Rabi oscillations occurring.

In the present work, we examine the temporal dynamics of a photon Bose-Einstein condensate (BEC) in the strongly fluctuating regime, where a steady drivendissipative state is induced by a balance of continuous dye pumping and cavity losses. We observe distinct temporal oscillations of the photon number correlations g (2) ( τ ),

even though the spectral photon distribution is, within experimental uncertainties, indistinguishable from predictions for thermodynamic equilibrium (Bose-Einstein distribution). Temporal oscillations of density fluctuations have been observed, e.g. in oscillatory relaxation dynamics in lasers [25] and in driven-dissipative atomic Bose-Einstein condensates confined in high-finesse optical cavities [26]. Damped, oscillatory displacement dynamics have been observed in colloids suspended in liquids upon driving out of equilibrium [27], see also theoretical work proposing corresponding experiments with laser-driven quantum dots [28]. In all those systems, however, the (temporally averaged) spectrum clearly differs from the thermodynamic equilibrium distribution.

For the correlation dynamics of the photon BEC, we find quantitative agreement with a theoretical analysis in terms of nonlinear rate equations. We are thus able to trace back the origin of the correlation oscillations to the coupling of the dye reservoir and the condensate photons, combined with time-reversal symmetry breaking in the (weakly) driven-dissipative state. The nonlinearity of the coupling implies that the dynamics of g (2) ( τ ) will depend on the steady state of the system.

The paper is organized as follows. Section II provides some details of the experimental setup and mode of operation. In Section III we derive the rate equations for the nonequilibrium dynamics of the average numbers of photons and of dye-molecule excitations as well as for the autocorrelation functions of these quantities. In Section IV we present and discuss the experimental results, along with the comparison to the theory. We conclude in Section V with an outlook to further studies.

## II. EXPERIMENTAL SETUP AND MODE OF OPERATION

Our experimental setup for trapping a two-dimensional photon gas in a dye-filled optical microcavity is similar to that described in earlier work [7, 29], as shown in Fig. 1. The cavity is composed of two highly reflecting mirrors (reflectivity &gt; 99 998%) of 1 m curvature radius spaced . by ≈ 2 µ mdistance, and is filled with rhodamine dye dissolved in ethylene glycol (1 mmol/L). Due to the small mirror spacing, the cavity has a longitudinal mode spacing comparable to the emission width of the dye. In this regime, we observe that only photons of a fixed longitudinal mode number q populate the cavity, with q = 7 in our setup. This imposes an upper limit on the optical wavelength and a restriction of energies to a minimum cutoff of glyph[planckover2pi1] ω c ≈ 2 1 . eV for photons in the cavity, where ω c is the cutoff frequency of the cavity with transverse momentum k r = 0. The optical dispersion becomes quadratic and the mirror curvature induces a trapping potential for the photons (see Fig. 1). One can show that the system is formally equivalent to a harmonically trapped two-dimensional gas of massive particles for which - other than for (three-dimensional) blackbody

Figure 1: The photon flow through the condensed photon gas in the dye-filled microcavity. The dye molecules, which are excited with a pump rate Γ ↑ constitute a heat and particle reservoir for the condensed photons. The photon dispersion in the cavity in terms of the transverse wave vector is hν = √ ( hν c ) 2 +( glyph[planckover2pi1] ck r ) 2 , with ν c = ω / π c 2 the cavity cutoff frequency and c the light velocity in the dye medium. The photon loss rate of the condensate through the cavity mirrors is represented by κ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00883v2-figure-f0762606-f441-47c9-8b21-96307ba0cfde.jpg)

radiation - Bose-Einstein condensation is possible for a thermalized ensemble [7]. To initially populate the cavity and compensate for losses from e.g., mirror transmission, the dye is pumped to a stationary state with an external laser beam.

To record the statistics as well as the number correlations of the photon condensate, the microcavity emission is directed through an optical telescope and a mode filter consisting of a 10 6 mm diameter iris to separate the . condensate mode from the higher-mode (thermal cloud) photons. Even though the divergence of the higher-order modes is larger than that of the ground mode, they contribute residually to the intensity near the optical axis in the far field, and a fraction of photons in these modes is consequently transmitted through the momentum filter. This imperfection of the mode filtering reduces the magnitude of the measured (normalized) intensity correlations. As a result, the measured correlations are smaller than the calculated values for the respective system parameters. We do not expect the temporal shape of the correlation signal to be affected by the imperfect mode filtering.

The light passing through the momentum filter is then sent through a polarizer to remove the polarization degeneracy and then imaged onto a fast (5 GHz bandwidth) photomultiplier. The electronic output signal of the photomultiplier is analyzed employing a 3 5 GHz bandwidth . oscilloscope. To suppress the influence of electronic noise of the high-bandwidth electronic analysis system (which is mainly attributed to the oscilloscope's analog-digital converters), the photomultiplier output is simultaneously recorded by two oscilloscope channels, with the crosscorrelation used for further data analysis of the fluctuations. Calibration of the photomultiplier signal is performed via the measured spectra, which relate the photon number in the condensate peak to the known

critical photon number in the thermal photon cloud of N c = π / 2 3( k B T/ glyph[planckover2pi1] Ω) 2 = 80660 for the used experimental parameters ( T = 300 K, trap frequency Ω = 2 π · 40 GHz).

The dephasing time τ ϕ of dye excitations (given by the collision time of solvent molecules with the dye on a scale of 10 -14 s is much shorter than the photon lifetime in the cavity ( τ p ≈ 0 5 . ns, determined by cavity losses), τ ϕ glyph[lessmuch] τ p . Furthermore, coherence between dye excitations and photons cannot be established, that is, our experiment operates in the weakly coupled regime of matter and light, and the trapped particles are photons, not polaritons [23, 24]). The photon lifetime, in turn, is much shorter than the nonradiative decay time of dye excitations, τ p glyph[lessmuch] τ nrad (with τ nrad being in the order of 50 ns [30]. This means that the sum of the number of cavity photons (including the thermal mode photons) and the dye excitation number can be considered conserved on the time scale of the photon lifetime τ p , where the photon number alone is fixed on average only. To good approximation, also the sum of the number of condensate mode photons and dye excitations is conserved, as exchange between thermal cloud photons and condensate mode photons only occurs via the bath and the dye excitation number exceeds the (thermal mode) photon number by many orders of magnitude [12] For sufficiently fast energy exchange with the dye reservoir, i.e., if several absorption and re-emission cycles occur before a photon is lost, the photons reach a thermal spectral distribution at the rovibrational temperature of the dye, which is at room temperature (300 K). Fig. 2 shows measurements of the optical spectrum for different average condensate photon numbers 〈 n 〉 ∞ (different pump rates). The visible good agreement of the the experimental data (dots) with the predictions of equilibrium theory (solid lines) confirms that in our experiments, a thermal distribution is achieved to good approximation, despite continuous pumping and losses. Above a critical photon number, the thermal photon gas eventually forms a BEC, signaled by the condensate peak in the spectral distribution at the position of the low-frequency cavity cutoff on top of a broad thermal photon cloud [7].

By varying the photon number 〈 n 〉 ∞ with respect to the number of dye reservoir molecules M within the cavity volume, the photon number statistics can be continuously tuned from a small relative reservoir size with Poissonian statistics and g (2) (0) ≈ 1 to a strongly fluctuating state of small relative photon number with Bose-Einstein statistics and g (2) (0) ≈ 2 [21].

## III. THEORETICAL DESCRIPTION

## A. Model and master equation

The thermalization dynamics of the dye-filled microcavity has been studied [19, 31, 32] employing a FrankCondon model for the dye molecules: The two-level sys-

Figure 2: Optical spectra (arbitrary units) as a function of wavelength λ , taken from cavity emission through the mirrors at the cutoff wavelength of λ c = 571 3 nm, where the . experiment was conducted. The solid lines are fits of the Bose-Einstein distribution at 300K (broadened by the experimental resolution) to the experimental data, with the steadystate average condensate photon number 〈 n 〉 ∞ as the only adjustable parameter. The fitted values of 〈 n 〉 ∞ are shown in the legend. The spectra at different 〈 n 〉 ∞ are vertically shifted for clarity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00883v2-figure-8ad6a563-3859-4cce-a592-111375e53e96.jpg)

tem of the electronic ground and excited states, between which the optical transitions occur, is coupled to a molecular vibrational degree of freedom (phonon) because the rest position of the ionic molecular oscillator depends on whether the molecule is in its electronic ground or excited state. The total number of dye molecules is M = M ↑ + M ↓ , where M ↑ and M ↓ are the number of molecules in the electronically excited or ground state, respectively. The Hamiltonian for the dye-filled microcavity reads thus [33]

$$\stackrel { \sim } { \text{the} } \quad & \quad \cdot \quad \cdot \quad \cdot \\ \text{the} \quad & \text{by} \quad & \quad H = \sum _ { k } \omega _ { k } a _ { k } ^ { \dagger } a _ { k } + \sum _ { m = 1 } ^ { M } \left [ \frac { \delta } { 2 } \sigma _ { m } ^ { z } + \Omega b _ { m } ^ { \dagger } b _ { m } \\ \text{to} \quad & \quad + \Omega \sqrt { S } \, \sigma _ { m } ^ { z } \left ( b _ { m } + b _ { m } ^ { \dagger } \right ) + g \sum _ { k } \left ( a _ { k } \sigma _ { m } ^ { + } + a _ { k } ^ { \dagger } \sigma _ { m } ^ { - } \right ) \right ], \\ \text{in} \quad & \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdots$$

in units such that glyph[planckover2pi1] = 1. The cavity-photon modes with transverse dispersion ω k are represented by the bosonic operators a k , a † k , the vibronic states of dye molecule m with oscillator frequency Ω by b m , b † m , and the electronic two-level system of molecule m by the Pauli matrix σ z m and and raising/lowering operators σ ± m , with the electronic transition frequency δ . The Frank-Condon electron-phonon coupling is parametrized by S , where the phonon position operator is ˆ x ∝ ( b m + b † m ) . The last term in Eq. (1) describes photon emission or absorption with the optical transition matrix element g , the smallest energy scale in the system. Since we consider photon gases which have already reached a stationary, thermal

distribution, and since in the experiment the transverse cavity ground-state mode (condensate mode) is singled out, the analysis may be restricted to the photon correlations in the condensate, i.e., we will collapse the sum over cavity modes in Eq. (1) to k = 0, where ω c ≡ ω 0 is the cavity-cutoff frequency (see also Fig. 1). The cutoff is chosen such that the cavity detuning is ∆ = ω c -δ &lt; 0. From here on, we will drop the cavity-mode subscript on the photon operators and write a ≡ a 0 .

The molecular part of the Hamiltonian can be diagonalized by a polaron transformation [33]. This leads to an effective, nonlinear electron-photon coupling, mediated by the phonon excitations of the dye. Due to fast collisions of the dye molecules with solvent molecules, the phonon excitations may be considered to be in thermal equilibrium at ambient temperature. After treating the phonon excitations as a Markovian thermal bath, this coupling can be parametrized by a coherent part g β [34], and incoherent, phonon-assisted couplings B abs for photon absorption and B em for photon emission (see Fig. 1). For the density matrix ρ , one obtains in this way the master equation

$$\dot { \rho } = \mathfrak { i } \left [ \rho, H _ { 0 } \right ] + \kappa \mathcal { L } [ a ] \rho + \sum _ { m = 1 } ^ { M } \mathcal { L } _ { m } \rho, \quad \ \ ( 2 ) \quad \exp \nolimits _ { \text{diat} }$$

where H 0 = ∆ a a † + ∑ M m =1 g β ( a σ † -m + aσ + m ) is the Hamiltonian generating the coherent part of the evolution in the rotated frame. The parameter κ in Eq. (2) describes the cavity loss to the environment, where the Lindblad operator acting on the density matrix ρ is defined as L [ X ρ ] = 1 2 ( [ Xρ,X † ] + [ X,ρX † ] ) . The molecule-induced superoperator L m is given by

$$\mathcal { L } _ { m } = \Gamma _ { \uparrow } \mathcal { L } [ \sigma _ { m } ^ { + } ] + \Gamma _ { \downarrow } \mathcal { L } [ \sigma _ { m } ^ { - } ] + B _ { \text{abs} } \mathcal { L } [ a \sigma _ { m } ^ { + } ] + B _ { \text{em} } \mathcal { L } [ a ^ { \dagger } \sigma _ { m } ^ { - } ]. \quad \text{in the} \\ ( 3 ) \quad \text{limit},$$

The four terms in L m describe, in order of appearance, pumping by an external laser source (Γ ), nonradiative ↑ decay of dye excitations (Γ ↓ glyph[lessmuch] Γ ), and photon absorp-↑ tion ( B abs ) and emission ( B em ) by the dye molecules, respectively.

Since the experiment operates in the regime g β /g glyph[lessmuch] 1 [34], and the detuning is very large compared to the renormalized coherent coupling, | ∆ | glyph[greatermuch] g β , the cavity mode effectively couples to the molecules only incoherently via the Lindblad terms proportional to B abs or B em . Therefore, the contribution of H 0 to the master equation (2) can be neglected for the present setup [19]. In addition, one should note that B em &gt; B abs because of the red detuning of the cavity cutoff ω c with respect to the electronic dye excitation energy δ .

## B. Average particle numbers

The coupled rate equations for the average condensate mode photon number 〈 n 〉 and the number of dye glyph[negationslash]

molecules in excited states 〈 M ↑ 〉 can now be derived from the master equation, where 〈·〉 = Tr[ ρ t ( ) ] · denotes the thermal and quantum mechanical average. Inserting ˙ ρ from Eqs. (2) and (3), using cyclic permutation under the trace, and σ ± m m σ ∓ = (1 ± σ z m ) / 2, leads to operator products of M ↑ = ∑ m (1 + σ z m ) / 2 and n = a a † . Again because of the fast dye-solvent collisions, coherent propagation of excitations of different dye molecules ( m ′ = m ) is negligible. This means that the sum over a large number of molecules amounts to an average, 〈 M ↑ 〉 = ∑ m (1 + 〈 σ z m 〉 ) / 2 = M (1 + 〈 σ z m 〉 ) / 2, and expectation values of higher-order operator products factorize, 〈 nM ↑ 〉 ≈ 〈 n 〉〈 M ↑ 〉 . In this way, one obtains the nonlinear, coupled rate equations

$$\underset { \substack { \ell \\ \ell } } { \ell } \quad \begin{array} { c } \ell & \ell \\ \frac { d } { d t } \langle n \rangle = - \kappa \langle n \rangle - B _ { \text{abs} } \left ( M - \langle M _ { \uparrow } \rangle \right ) \langle n \rangle \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array} \quad \begin{array} { c } \ell \\ \ell \\ \ell } \end{array}$$

$$\mathfrak { h } _ { \Phi } ^ { \dagger } \left ( M - \left \langle M _ { \uparrow } \right \rangle & = \Gamma _ { \uparrow } \left ( M - \left \langle M _ { \uparrow } \right \rangle \right ) - \Gamma _ { \downarrow } \left \langle M _ { \uparrow } \right \rangle \\ \mathfrak { h } _ { \Phi } ^ { \dagger } \left ( M - \left \langle M _ { \uparrow } \right \rangle \right ) - \Gamma _ { \downarrow } \left \langle M _ { \uparrow } \right \rangle \\ \mathfrak { h } _ { \Phi } ^ { \dagger } \left ( M - \left \langle M _ { \uparrow } \right \rangle \right ) \left \langle n \right \rangle - B _ { \text{em} } \left ( \left \langle n \right \rangle + 1 \right ) \left \langle M _ { \uparrow } \right \rangle. \\ & \quad \right. \quad \dots \quad \dots \quad \dots \quad \dots.$$

These agree, in fact, with the semiclassical rate equations expected phenomenologically from pumping and nonradiative decay of molecule excitations a well as stimulated and spontaneous photon emission into the cavity. Alternatively, one may solve the untruncated rate equations for 〈 n 〉 and 〈 M ↑ 〉 together with three equations for the second moments [31]. This is discussed in detail in the Appendix. For large M , both solution methods give the same results for 〈 n 〉 and 〈 M ↑ 〉 in the long-time limit.

$$\langle n \rangle _ { \infty } = \frac { M \left ( B _ { \text{em} } \Gamma _ { \uparrow } - B _ { \$$

To calculate the steady-state second-order photon correlations in the next section, it will be necessary to have the average numbers of photons and excited molecules in the steady state which is reached in the long-time limit, 〈 n 〉 ∞ and 〈 M ↑ 〉 ∞ . This amounts to setting the time derivatives in Eqs. (4) to zero, and one obtains for large molecule number M glyph[greatermuch] 1,

$$\langle M _ { \uparrow } \rangle _ { \infty } = \frac { M B _ { \text{abs} } + \kappa } { B _ { \text$$

In our experiments, the pump rate strongly exceeds the nonradiative decay, Γ ↑ glyph[greatermuch] Γ , ↓ and B em + B abs ≈ B em . The ratio of emission and absorption is given by B em /B abs = exp( -glyph[planckover2pi1] ∆ /k B T ), where T is the phonon temperature (Kennard-Stepanov relation). With these simplifications, the steady-state photon number becomes approximately

$$\langle n \rangle _ { \infty } \approx \frac { M \Gamma _ { \uparrow } } { \kappa }. \quad \quad \quad$$

This expression is useful for converting 〈 n 〉 ∞ , which is measured in the experiments, into the pump parameter Γ ↑ of the theoretical model and vice versa. When comparing to experimental data, however, a full numerical solution for the steady state of Eqs. (4) and (25) is used.

## C. Second-order correlation function

The time-dependent photon density-density or secondorder correlation function measured in the experiment is defined as

$$g ^ { ( 2 ) } ( \tau ) = \frac { \langle n ( t + \tau ) n ( t ) \rangle } { \langle n ( t ) \$$

where ˆ is the total Liouvillian superoperator belonging L to the master equation (2), ρ ∞ = lim t →∞ ρ t ( ) denotes the steady-state density matrix, and we define ˜ ρ ∞ := a aρ † ∞ . Note in passing that, for the normal-ordered second-order correlation function, one would need to set ˜ ρ ∞ = aρ ∞ † a . Defining also an effective average X := Tr X e ˆ Lτ ˜ ρ ∞ , one has g (2) ( τ ) = n/ n 〈 〉 2 ∞ . Formally, n and M ↑ obey almost the same definitions as 〈 n 〉 and 〈 M ↑ 〉 , however with ρ t ( ) replaced by e ˆ Lτ ˜ ρ ∞ and Tr ˜ ρ ∞ = 〈 n 〉 ∞ . Thus, one finds equations of motion analogous to Eqs. (4a) and (4b),

$$\frac { \text{d} } { \text{d} \tau } \overline { n } & = - \kappa \overline { n } - B _ { \text{$$

$$\frac { \cdot - \infty \cdot \cdot \cdot \cdot - \infty - 1 } { d \tau } \overline { M _ { \uparrow } } & = \Gamma _ { \uparrow } \left ( \langle n \rangle _ { \infty } M - \overline { M _ { \uparrow } } \right ) - \Gamma _ { \downarrow } \overline { M _ { \uparrow } } \\ & + B _ { \text{abs} } \overline { ( M - M _ { \uparrow } ) n } - B _ { \text{em} } \overline { ( n + 1 ) M _ { \uparrow } }. \quad \ ( 9 b ) \quad \text{The } \cdot \\ & \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The cluster expansions of averages of the higher-order operator products nM ↑ are truncated under the same conditions as discussed in the previous section (incoherent propagation of excitations of different molecules). This amounts to letting nM ↑ = 〈 n 〉 ∞ M ↑ + 〈 n 〉 ∞ 〈 nM ↑ 〉 ∞ + 〈 M ↑ 〉 ∞ n - 〈 〉 2 n 2 ∞ 〈 M ↑ 〉 ∞ . By computing the steady-state density matrix numerically exactly for different molecule numbers of order M ∼ 10 , 2 we have checked that the truncation is well justified already for intermediate system sizes, i.e. above threshold [31]. The experiment is operated far inside this range of validity. Above threshold, we also have Γ ↑ glyph[greatermuch] Γ , ↓ such that we may safely neglect the non-radiative molecule decay rate Γ ↓ for the evaluations in Section IV.

In terms of the deviations of the second-order correlation functions from their relaxed values (attained for τ →∞ ),

$$g = \begin{pmatrix} \Delta g _ { n } ^ { ( 2 ) } \\ \Delta g _ { n, \, M _ { \uparrow } } ^ { ( 2 ) } \end{pmatrix} = \begin{pmatrix} \overline { n } - \langle n \rangle _ { \infty } ^ { 2 } \\ \overline { M } _ { \uparrow } - \langle n \rangle _ { \infty } \langle M _ { \uparrow } \rangle _ { \infty } \end{pmatrix}, \quad ( 1 0 )$$

Eqs. (9) then become

$$\text{eq.} & \text{qs.} \ ( 9 ) \text{ then become} \\ \text{ad-t is} \\ & + B _ { \text{em} } \left ( \langle n \rangle _ { \infty } + 1 \right ) \langle M _ { \uparrow } \rangle _ { \infty } \right ] - \kappa \Delta g _ { n } ^ { ( 2 ) } \\ & - B _ { \text{abs} } \left [ \left ( M - \langle M _ { \uparrow } \rangle _ { \infty } \right ) \Delta g _ { n } ^ { ( 2 ) } - \langle n \rangle _ { \infty } \Delta g _ { n, M _ { \uparrow } } ^ { ( 2 ) } \right ] \\ & + B _ { \text{em} } \left [ \left ( \langle n \rangle _ { \infty } + 1 \right ) \Delta g _ { n, M _ { \uparrow } } ^ { ( 2 ) } + \langle M _ { \uparrow } \rangle _ { \infty } \Delta g _ { n } ^ { ( 2 ) } \right ], \\ \frac { d } { } \Delta a ^ { ( 2 ) }. \tau \ & = \langle n \rangle _ { \infty } \left [ \Gamma _ { \ast } \left ( M - \langle M _ { \ast } \rangle _ { \infty } \right ) - \Gamma _ { \ast } \langle M _ { \ast } \rangle _ { \infty }$$

$$& \cdot, \ \vee \nwarrow & & ( 1 1 a ) \\ & \frac { d } { d \tau } \Delta g _ { n, M _ { \uparrow } } ^ { ( 2 ) } & = \langle n \rangle _ { \infty } \left [ \Gamma _ { \uparrow } \left ( M - \langle M _ { \uparrow } \rangle _ { \infty } \right ) - \Gamma _ { \downarrow } \langle M _ { \uparrow } \rangle _ { \infty } \\ \intertext { n } _ { \infty } & \coloneqq \\ & \frac { + } { \infty } \coloneqq \\ & \frac { - } { \infty } \coloneqq \\ \text{ded} \quad \text{to set} \quad & \frac { + } { \infty } \right ) \Delta g _ { n } ^ { ( 2 ) } - \langle n \rangle _ { \infty } \Delta g _ { n, M _ { \uparrow } } ^ { ( 2 ) } \right ] \\ & \frac { + } { \infty } \coloneqq \\ & \frac { - } { \nmid } \end{matrix} & \frac { - } { \infty } \left [ \langle \langle n \rangle _ { \infty } + 1 \Delta g _ { n, M _ { \uparrow } } ^ { ( 2 ) } + \langle M _ { \uparrow } \rangle _ { \infty } \Delta g _ { n } ^ { ( 2 ) } \right ]. \\ & \langle M _ { \uparrow } \rangle, \quad & \text{Using the steady-state solution of Eqs. } ( 4 ), \, \text{one eventu}$$

Using the steady-state solution of Eqs. (4), one eventually finds a system of two coupled linear equations,

$$\partial _ { \tau } g & = \begin{pmatrix} - \kappa - \tilde { \Gamma } _ { M } & \tilde { \Gamma } _ { n } \\ \tilde { \Gamma } _ { M } & - ( \Gamma _ { \uparrow } + \Gamma _ { \downarrow } ) - \tilde { \Gamma } _ { n } \end{pmatrix} g. \quad \text{ (12)}$$

The matrix elements are given by

$$\tilde { \Gamma } _ { M } & = B _ { \text{abs} } ( M - \langle M _ { \uparrow } \rangle _ { \infty } ) - B _ { \text{em} } \langle M _ { \uparrow } \rangle _ { \infty }, \\ \tilde { \Gamma } _ { n } & = B _ { \text{abs} } \langle n \rangle _ { \infty } + B _ { \text{em} } ( \langle n \rangle _ { \infty } + 1 ). \\ \tau _ { \infty } \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

The coupling constant ˜ Γ M is composed of an absorption term proportional to the number of ground-state molecules in the steady state ( M -〈 M ↑ 〉 ∞ ), and a corresponding emission term with the number of excited molecules. The coupling constant Γ ˜ n is given by an absorption term, and the terms corresponding to stimulated and spontaneous emission.

The result of Eq. (12) is equivalent to what one would obtain from linearizing Eqs. (4) around the steady state and then applying the regression theorem.

## IV. RESULTS

The coupling matrix in Eq. (12) is non-Hermitian because of time-reversal symmetry breaking in the drivendissipative system. As a result, its eigenvalues are found to be complex, λ = λ ′ ± i λ ′′ , where λ , λ ′ ′′ ∈ R . In terms of the quantities defined in Eqs. (13), the eigenvalues are given by

$$\lambda = - \gamma \pm \sqrt { \gamma ^ { 2 } - \omega _ { 0 } ^ { 2 } },$$

where

$$2 \gamma = \kappa + \Gamma _ { \uparrow } + \Gamma _ { \downarrow } + \tilde { \Gamma } _ { M } + \tilde { \Gamma } _ { n }, \quad \ \ ( 1 5 a )$$

$$\omega _ { 0 } ^ { 2 } = \kappa \tilde { \Gamma } _ { n } + ( \kappa + \tilde { \Gamma } _ { M } ) ( \Gamma _ { \uparrow } + \Gamma _ { \downarrow } ). \quad ( 1 5 b )$$

Figure 3: Typical results for the time dependence of the second-order correlation function for two different photon numbers in the condensate, as indicated. The solid lines are fits to the model function, Eq. (19). The envelope of the observed correlation signal g (2) ( τ ) temporally decays with a typical relaxation time of about 4 ns.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00883v2-figure-7aa6d974-33bf-47c0-ab33-2514c741b606.jpg)

For Γ ↓ = 0, B abs = αB em and to leading order in 1 /M , this reduces to more insightful approximate expressions. In this case, we have

$$\langle n \rangle _ { \infty } \approx M \frac { \Gamma _ { \uparrow } } { ( 1 + \alpha ) \kappa }, \quad \quad ( 1 6 \tt a ) \quad \underset { \tau _ { b } } { \sum } \sum _ { \tau _ { b } } ^ { \ m a }$$

$$\langle M _ { \uparrow } \rangle _ { \infty } \approx \frac { M \alpha B _ { \text{em} } + \kappa } { ( 1 + \alpha ) B _ { \text{em} } }. \quad \quad ( 1 6 b ) \quad \text{$\text{$Eq.}$} \\ \text{$\text{no}$}$$

This gives Γ ˜ M = -κ and ˜ Γ n = M Γ ↑ B em /κ , from which we find

$$\gamma = \frac { \Gamma _ { \uparrow } } { 2 } \left ( 1 + \frac { M B _ { \text{em} } } { \kappa } \right ) \approx \frac { M \Gamma _ { \uparrow } B _ { \text{em} } } { 2 \kappa }, \quad ( 1 7 \text{a} ) \quad \text{relat} \\ \text{$\quad$func$} \quad \text{func$} \quad \text{in $Fi$}$$

$$\omega _ { 0 } ^ { 2 } & = M \Gamma _ { \uparrow } B _ { \text{em} }. & ( 1 7 b ) & \quad \overset { \mu \, \Gamma \, \iota } { \text{the c} }$$

Hence, the eigenvalues are given by

$$\lambda = - \frac { M \Gamma _ { \uparrow } B _ { \text{em} } } { 2 \kappa } \pm \frac { 1 } { 2 } \sqrt { \left ( \frac { M \Gamma _ { \uparrow } B _ { \text{em} } } { \kappa } \right ) ^ { 2 } - 4 M \Gamma _ { \uparrow } B _ { \text{em} } }. \quad \begin{array} { c } \text{$\text{$seen$} \\ $numb$} \\ $the t l \\ $ \text{perim$} \\ $ R \dots \, \end{array} }$$

Without drive and dissipation (Γ ↑ = κ = 0), we have a two-component system which shows a single relaxation time, i.e. the system possesses one zero eigenvalue and a purely real eigenvalue λ = -( ˜ Γ M + ˜ Γ ). Correspondingly, n the open-system character here is necessary to achieve an imaginary part for the eigenvalues and hence an oscillating photon-number correlation.

Note that this approximation leading to Eq. (18) is not employed in the analysis of the experimental data. It rather serves to illustrate how the openness of the system influences the eigenvalues through the steady-state occupations by making the dependence on Γ ↑ and κ explicit.

For the second-order correlation function, one hence finds a solution of the form

$$g _ { n } ^ { ( 2 ) } ( \tau ) = 1 + e ^ { \lambda ^ { \prime } \tau } \left [ c _ { 1 } \cos ( \lambda ^ { \prime \prime } \tau ) + c _ { 2 } \sin ( \lambda ^ { \prime \prime } \tau ) \right ], \quad ( 1 9 )$$

with the real part λ ′ &lt; 0. The initial values for the dynamics of the second-order correlation functions are found from the steady-state solutions for the second moments, 〈 n 2 〉 ∞ and 〈 nM ↑ 〉 ∞ . Typical experimental data for the temporal variation of the second-order coherence function are given in Fig. 3 for average photon numbers 〈 n 〉 ∞ ≈ 4620 and 〈 n 〉 ∞ ≈ 17100, respectively, showing damped, oscillatory behavior, as expected from the theoretical analysis. Within the studied experimental parameters, the oscillations persist as the system is tuned between grand canonical and canonical statistical regimes. As can be seen in the two examples of Fig. 3, the magnitude of photon number fluctuations (relative to the average photon number) decreases at higher average photon numbers corresponding to a smaller relative size of the dye reservoir. This behavior has been studied in more detail in earlier work [21]. For the present system, the theoretical values for g (2) (0) of the condensate mode at 〈 n 〉 ∞ ≈ 4620 and 〈 n 〉 ∞ ≈ 17100 are ≈ 2 and ≈ 1 3, . respectively. As described above, the smaller g (2) (0) values observed in the experimental system are attributed to imperfect mode filtering, which tends to reduce the magnitude of the correlation signal towards that of an uncorrelated sample, where g (2) ( τ ) = 1 for all times τ . The solid lines in Fig. 3 are fits of the function from Eq. (19) to the experimental data, with the imperfect mode filtering resulting in a lowering of the coefficients c 1 , c 2 with respect to the theory values. The experimental values of the second-order relaxation time τ c = 1 / λ | ′ | and the oscillation frequency ω (2) = λ ′′ of the g (2) correlations are determined by fitting the theoretical model function (Eq. (19)) to the experimental data, as depicted in Fig. 3. In this way, we have recorded the variation of the oscillation frequency ω (2) upon the change of the average photon number 〈 n 〉 ∞ , as shown in Fig. 4 (dots). We observe an increase of the oscillation frequency of the second-order coherence function with the average photon number. The solid line in Fig. 4 is obtained using a fit of the theoretical eigenvalues λ ′ , λ ′′ of Eq. (12) to the experimental data, where the model parameters κ , Γ , and ↑ B em were used as fit parameters, and the nonradiative decay rate Γ ↓ was set to zero.

The experimental data are fitted to good precision for all different 〈 n 〉 ∞ by three parameters which are consistent with experimentally estimated values. We interpret this (as well as the comparison shown in Fig. 3)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00883v2-figure-35f7795b-c7b5-4109-89df-a8d8e22acccb.jpg)

×

Figure 4: Oscillation frequency ω (2) of the second order correlation function g (2) ( τ ) as a function of the average photon number 〈 n 〉 ∞ in the condensate, as measured in the experiment (dots) and predicted by a nonlinear rate equation model (solid line). See the main text for details. The parameter values are M = 5 17 . · 10 , 9 κ = 2 33 GHz, . B em = 2 50 . · 10 -5 GHz, B em /B abs ≈ 57 (Kennard-Stepanov relation corresponding to the cutoff wavelength of λ c = 571 3 nm), Γ . ↓ = 0.

as evidence that the origin of the g (2) oscillations can be traced back to the effects incorporated in our rate equation model, namely time-reversal symmetry breaking due to nonequilibrium pumping and dissipation, and the coupling between the subsystems of dye-molecule excitations and cavity photons.

Remarkably, despite the clear nonequilibrium signatures in g (2) ( τ ), a spectral distribution which, within experimental accuracies, is indistinguishable from the thermal equilibrium Bose-Einstein distribution is attained (see Fig. 2). The physical origin of this seemingly contradictory behavior is that there is a (stationary) net flow of photons from the dye reservoir to the resonator and out to the environment, in such a way that the average photon number is constant and the photon gas reaches a thermal spectral distribution due to the thermal contact with the dye. Interestingly, our calculations show that the nonvanishing zero-delay second order correlations g (2) (0) &gt; 1, from the grand canonical nature of the photon condensate, are (in the investigated limits) not affected by the nonequilibrium character of the system. However, the residual photon loss and pump do affect the temporal decay of the correlation function g (2) ( τ ). This photon flow causes additional structure in the intensity correlations, as observed in the present work. A related behavior is known in nanoelectronic systems, e.g. current carrying, metallic nanowires, where the electron energy distribution is equilibrated from fast electron-electron collisions in the wire [35].

## V. CONCLUSIONS

To conclude, we observed an oscillatory behavior of grand canonical Bose-Einstein condensates by studying the second-order coherence of the emission of a dye microcavity. Its origin is traced back to the remnant drivendissipative character of the light condensate. Our results show that even when the energy distribution of particles to good accuracy follows the predictions for thermal equilibrium, fluctuation dynamics depend sensitively on the openness of the system. We note that a related behavior can be observed in the hot-electron regime of electronic quantum wires at large bias voltage, where nonthermal noise, albeit not oscillatory, is induced by the current. Due to fast electron-electron collisions in the wire, it coexists with an equilibrium (Fermi-Dirac) distribution of the electron energy. Here we observe this phenomenon in a photon system for the first time.

The g (2) ( τ ) oscillations observed in our photon condensates are reminiscent of relaxation oscillations in lasers. However, there are important differences. First, a laser is in a state far from equilibrium with nonthermal spectral distribution. In contrast, our system is operated in a near-equilibrium state with a thermal Bose-Einstein distribution, but clearly nonequilibrium dynamics of the photon number correlations. Secondly, in a steady-state laser g (2) ( τ ) = 1 for all times. This means that oscillations of g (2) ( τ ) do not occur unless the system is perturbed from outside, such as in the initial increase of the cavity photon number when a laser is turned on. For our system, when it is not operated with a very small relative size of the reservoir, the steady state is characterized by g (2) (0) &gt; 1. These fluctuations of the grand canonical system are responsible for the excitation of the oscillations.

Our findings open up new avenues for further investigations of the open-system dynamics of grand canonical photon condensates. In the future, the experiments may be extended to further study the regime with stronger dissipation and drive. In addition, the second-order correlations can be used as a tool to sensitively characterize the system parameters. In a lattice with several coupled grand canonical photon condensates, a variety of new dynamical phases may be expected.

Note added in proof. Recently, a report on related experiments with a dye-microcavity subject to short pulsed picosecond laser pump irradiation appeared [36].

## ACKNOWLEDGMENTS

We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) within the Cooperative Research Center SFB/TR 185 (277625399) and the Cluster of Excellence ML4Q (390534769), from the European Union within the ERC project INPEC and the Quantum Flagship project PhoQuS, and from the DLR within project BESQ.

## APPENDIX: TRUNCATION OF HIERARCHY

In this Appendix, we present details on the hierarchy of the equations of motion for the expectation values of successively increasing order. These will be calculated from the master equation [19]

$$\dot { \rho } & = \frac { \kappa } { 2 } \mathcal { L } [ a ] \rho + \frac { 1 } { 2 } \sum _ { m = 1 } ^ { M } \left \{ \Gamma _ { \uparrow } \mathcal { L } [ \sigma _ { m } ^ { + } ] + \Gamma _ { \downarrow } \mathcal { L } [ \sigma _ { m } ^ { - } ] \\ & + B _ { \text{em} } \mathcal { L } [ a ^ { \dagger } \sigma _ { m } ^ { - } ] + B _ { \text{abs} } \mathcal { L } [ a \sigma _ { m } ^ { + } ] \right \} \rho.$$

Assuming that the molecules are all identical, one can replace the sum over the M molecules in the equation for the photon occupation 〈 a a † 〉 by a factor of M and find

$$\partial _ { t } \langle a ^ { \dagger } a \rangle = - \kappa \langle a ^ { \dagger } a \rangle - \frac { M } { 2 } B _ { \text{abs} } \langle a ^ { \dagger } a \left ( 1 - \sigma ^ { z } \right ) \rangle$$

$$\partial _ { t } \langle a ^ { \prime } a \rangle & = - \kappa \langle a ^ { \prime } a \rangle - \frac { \mu } { 2 } B _ { \text{abs} } \langle a ^ { \prime } a \left ( 1 - \sigma ^ { z } \right ) \rangle \\ & + \frac { M } { 2 } B _ { \text{em} } \langle a a ^ { \dagger } \left ( 1 + \sigma ^ { z } \right ) \rangle, \\ \partial _ { t } \langle \sigma ^ { z } \rangle & = \Gamma _ { \uparrow } \left ( 1 - \langle \sigma ^ { z } \rangle \right ) - \Gamma _ { \downarrow } \left ( 1 + \langle \sigma ^ { z } \rangle \right ) \\ & + B _ { \text{abs} } \langle a ^ { \dagger } a \left ( 1 - \sigma ^ { z } \right ) \rangle - B _ { \text{em} } \langle a a ^ { \dagger } \left ( 1 + \sigma ^ { z } \right ) \rangle, \\ \dots \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \alpha \left ( \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Lambda \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \beta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Depsilon \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Dots \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta$$

where we have dropped the molecule index m . Multiplying Eq. (21b) by M/ 2 and using M∂ σ t 〈 z 〉 / 2 = ∂ t 〈 M ↑ 〉 , 〈 σ + σ -〉 = 1 + 〈 σ z 〉 / 2 = 〈 m /M 〉 , 〈 σ -σ + 〉 = 1 〈 -σ z 〉 / 2 = 〈 M -m /M 〉 , and 〈 nM ↑ 〉 ≈ 〈 n 〉〈 M ↑ 〉 , we arrive at Eqs. (4). The latter approximation is excellent for large systems: while ∆ g (2) n, M ↑ (0) = 〈 M n ↑ 〉 ∞ -〈 〉 n ∞ 〈 M ↑ 〉 ∞ does not vanish in general, its influence on 〈 n 〉 ∞ , 〈 M ↑ 〉 ∞ turns out to be negligible.

Accordingly, the equations of motion of the next-order expectation values in Eqs. (21) are given by

$$\partial _ { t } ( a ^ { \dagger } a ^ { \dagger } a a ) & = - 2 \kappa ( a ^ { \dagger } a ^ { \dagger } a a ) - 2 M B _ { \text{abs} } ( a ^ { \dagger } a a \sigma ^ { - } \sigma ^ { + } ) + 2 M B _ { \text{em} } \left ( ( a ^ { \dagger } a a \sigma ^ { + } \sigma ^ { - } ) + 2 ( a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } ) \right ), \quad ( 2 2 a ) \\ \partial _ { t } ( a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } ) & = - \kappa ( a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } ) + \Gamma _ { \uparrow } ( a ^ { \dagger } a \sigma ^ { - } \sigma ^ { + } ) - \Gamma _ { \downarrow } ( a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } ) \\ & + B _ { \text{med} } ( a ^ { \dagger } a ^ { \dagger } a a \sigma ^ { - } \sigma ^ { + } ) - B _ { \dots } \left ( ( a ^ { \dagger } a ^ { \dagger } a a \sigma ^ { + } \sigma ^ { - } ) + 2 ( a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } ) \right )$$

-

-

〈

〉

-

〈

〉

$$& + B _ { \text{abs} } \langle a ^ { \dagger } a ^ { \dagger } a a \sigma ^ { - } \sigma ^ { + } \rangle - B _ { \text{em} } \left ( \langle a ^ { \dagger } a ^ { \dagger } a a \sigma ^ { + } \sigma ^ { - } \rangle + 2 \langle a ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } \rangle \right ) \\ & \quad \, / \lambda \, \int \, \mathfrak { m } \, \quad \, / \dots + \_ - \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \_ \end{cases}$$

(

M

1)

B

abs a aσ

†

+

σ

-

˜

σ

-

˜

σ

+

+(

M

1)

B

em aa σ

†

+

σ

-

˜

σ

+

˜

σ

-

,

$$\partial _ { t } ( \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } ) & = \Gamma _ { \uparrow } \left ( ( \sigma ^ { - } \sigma ^ { + } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } ) + ( \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { - } \tilde { \sigma } ^ { + } ) \right ) - 2 \Gamma _ { \downarrow } ( \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } ) \\ & + B _ { \text{abs} } \left ( ( \real ^ { \dagger } a \sigma ^ { - } \sigma ^ { + } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } ) + ( \real ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { - } \tilde { \sigma } ^ { + } ) \right ) - 2 B _ { \text{em} } ( \real ^ { \dagger } a \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } ).$$

The Pauli matrices ˜ σ ± describe any molecule that is not identical to σ ± . Under the assumption that the total density matrix of the molecules is an incoherent mixture of all states corresponding to an excitation number of M ↑ [31], which will be the case for the steady-state density matrix of the master equation (20), one can show that the expectation values of four Pauli matrices decompose as

$$M ( M - 1 ) \langle \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } \rangle & = \langle M _ { \uparrow } ^ { 2 } \rangle - \langle M _ { \uparrow } \rangle, \\ M ( M - 1 ) \langle \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { - } \tilde { \sigma } ^ { + } \rangle & = M \langle M _ { \uparrow } \rangle - \langle M _ { \uparrow } ^ { 2 } \rangle.$$

$$( 2 3 )$$

Then the truncation of Eqs. (22) can be performed rigorously by expanding the highest-order expectation values according to

$$& M ( M - 1 ) \langle a a ^ { \dagger } \sigma ^ { + } \sigma ^ { - } \tilde { \sigma } ^ { + } \tilde { \sigma } ^ { - } \rangle \\ & = \langle ( n + 1 ) ( M _ { \uparrow } ^ { 2 } - M _ { \uparrow } ) \rangle \\ & = \langle n M _ { \uparrow } ^ { 2 } \rangle - \langle n M _ { \uparrow } \rangle + \langle M _ { \uparrow } ^ { 2 } \rangle - \langle M _ { \uparrow } \rangle,$$

and using the relation

$$\langle n M _ { \uparrow } ^ { 2 } \rangle = 2 \langle M _ { \uparrow } \rangle \langle n M _ { \uparrow } \rangle + \langle n \rangle \langle M _ { \uparrow } ^ { 2 } \rangle - 2 \langle n \rangle \langle M _ { \uparrow } \rangle ^ { 2 }.$$

In this manner, after multiplying Eq. (22b) by a factor of M and Eq. (22c) by M M ( -1), one obtains [31]

- ∂ t 〈 n 2 〉 = ∂ t 〈 a aa a † † 〉 = ∂ t 〈 a a aa † † 〉 + ∂ t 〈 a a † 〉 = κ ( 〈 n 〉 -2 〈 n 2 〉 ) -B abs [ 2 〈 n 2 〉 ( M -〈 M ↑ 〉 ) + 4 〈 M ↑ 〉〈 n 〉 2 - 〈 〉〈 4 n nM ↑ 〉 - 〈 n M ( -M ↑ ) 〉 ] + B em [ 4 〈 n 〉〈 nM ↑ 〉 +2 〈 M ↑ 〉〈 n 2 〉 -4 〈 M ↑ 〉〈 n 〉 2 + (3 〈 n +1) M ↑ 〉 ] , (25a)
- = - 〈 κ nM ↑ 〉 +Γ ↑ 〈 n M ( -M ↑ ) 〉 -Γ ↓ 〈 nM ↑ 〉

∂ t 〈 nM ↑ 〉 = M∂ a aσ t 〈 † + σ -〉

- + B abs [ 2 〈 M ↑ 〉〈 n 〉 ( 〈 n 〉 - 〈 M ↑ 〉 ) -M n ( 〈 〉 + 〈 nM ↑ 〉 ) + 〈 n 〉〈 M 2 ↑ 〉 +2( 〈 M ↑ 〉 - 〈 n 〉 + 1 2 ) 〈 nM ↑ 〉 +( M -〈 M ↑ 〉 ) 〈 n 2 〉 ]

(25b)

- + B em [ 2 〈 M ↑ 〉〈 n 〉 ( 〈 n 〉 - 〈 M ↑ 〉 ) -〈 M ↑ 〉 +( 〈 n 〉 +1) 〈 M 2 ↑ 〉 +2( 〈 M ↑ 〉 - 〈 n 〉 -1) 〈 nM ↑ 〉 - 〈 M ↑ 〉〈 n 2 〉 ] ,

∂ t 〈 M 2 ↑ 〉 = M M ( -1) ∂ t 〈 σ + σ -˜ σ + ˜ σ -〉 + ∂ t 〈 M ↑ 〉

- + B abs ( (2 M -1) 〈 nM ↑ 〉 + M n 〈 〉 -4 〈 M ↑ 〉〈 nM ↑ 〉 -2 〈 n 〉〈 M 2 ↑ 〉 +4 〈 n 〉〈 M ↑ 〉 2 )
- = 2Γ ↑ ( M M 〈 ↑ 〉 + M -〈 M ↑ 〉 - 〈 M 2 ↑ 〉 ) +Γ ↓ ( 〈 M ↑ 〉 -2 〈 M 2 ↑ 〉 )

(25c)

- -B em ( 4 〈 M ↑ 〉〈 nM ↑ 〉 +2 〈 n 〉〈 M 2 ↑ 〉 -4 〈 n 〉〈 M ↑ 〉 2 +2 〈 M 2 ↑ 〉 -( 〈 n 〉 +1) 〈 M ↑ 〉 ) .

As mentioned in the main text, the steady-state solution of these equations is required for the initial values of the dynamics of the second-order correlation functions.

- [1] L. P. Pitaevski˘ ı and S. Stringari, Bose-Einstein Condensation , International Series of Monographs on Physics, Vol. 116 (Clarendon Press, Oxford, 2003).
- [18] Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Phys. Rev. Lett. 118 , 016602 (2017).
- [2] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80 , 885 (2008).
- [3] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269 , 198 (1995).
- [4] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75 , 3969 (1995).
- [5] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szyma´ nska, R. Andr´, J. L. Staehli, V. Savona, P. B. Lite tlewood, B. Deveaud, and S. Le Dang, Nature 443 , 409 (2006).
- [6] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Science 316 , 1007 (2007).
- [7] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Nature 468 , 545 (2010).
- [8] J. Marelic and R. A. Nyman, Phys. Rev. A 91 , 033813 (2015).
- [9] S. Greveling, K. L. Perrier, and D. van Oosten, Phys. Rev. A 98 , 013810 (2018).
- [10] I. Fujiwara, D. ter Haar, and H. Wergeland, J. Stat. Phys. 2 , 329 (1970).
- [11] R. M. Ziff, G. E. Uhlenbeck, and M. Kac, Phys. Rep. 32 , 169 (1977).
- [12] J. Klaers, J. Schmitt, T. Damm, F. Vewinger, and M. Weitz, Phys. Rev. Lett. 108 , 160403 (2012).
- [13] D. N. Sob'yanin, Phys. Rev. E 85 , 061120 (2012).
- [14] A. Zamora, L. M. Sieberer, K. Dunnett, S. Diehl, and M. H. Szyma´ska, Phys. Rev. X n 7 , 041006 (2017).
- [15] D. Vorberg, R. Ketzmerick, and A. Eckardt, Phys. Rev. A 97 , 063621 (2018).
- [16] H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. 82 , 1489 (2010).
- [17] T. Byrnes, N. Y. Kim, and Y. Yamamoto, Nature Physics 10 , 803 (2014).
- [19] P. Kirton and J. Keeling, Phys. Rev. Lett. 111 , 100404 (2013).
- [20] J. Schmitt, T. Damm, D. Dung, F. Vewinger, J. Klaers, and M. Weitz, Phys. Rev. A 92 , 011602(R) (2015).
- [21] J. Schmitt, T. Damm, D. Dung, F. Vewinger, J. Klaers, and M. Weitz, Phys. Rev. Lett. 112 , 030401 (2014).
- [22] J. R. Lakowicz, Principles of Fluorescence Spectroscopy , 3rd ed. (Springer, Boston, MA, 2006).
- [23] E. D. Angelis, F. D. Martini, and P. Mataloni, J. Opt. B 2 , 149 (2000).
- [24] H. Yokoyama and S. D. Brorson, J. Appl. Phys. 66 , 4801 (1989).
- [25] N. Takemura, J. Omachi, and M. Kuwata-Gonokami, Phys. Rev. A 85 , 053811 (2012).
- [26] F. Brennecke, R. Mottl, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Proc. Nat. Acad. Sci. USA 110 , 11763 (2013).
- [27] J. Berner, B. M¨ uller, J. R. Gomez-Solano, M. Kr¨ger, u and C. Bechinger, Nat. Commun. 9 , 999 (2018).
- [28] T. Moradi, M. B. Harouni, and M. H. Naderi, Sci. Rep. 8 , 12435 (2018).
- [29] J. Klaers, J. Schmitt, T. Damm, F. Vewinger, and M. Weitz, Appl. Phys. B 105 , 17 (2011).
- [30] J. Barroso, A. Costela, I. Garcia-Moreno, and R. Sastre, Chem. Phys. 238 , 257 (1998).
- [31] P. Kirton and J. Keeling, Phys. Rev. A 91 , 033826 (2015).
- [32] J. Keeling and P. Kirton, Phys. Rev. A 93 , 013829 (2016).
- [33] M. Marthaler, Y. Utsumi, D. S. Golubev, A. Shnirman, and G. Sch¨n, Phys. Rev. Lett. o 107 , 093901 (2011).
- [34] M. Radonji´ c, W. Kopylov, A. Balaˇ, z and A. Pelster, New J. Phys. 20 , 055014 (2018).
- [35] A. H. Steinbach, J. M. Martinis, and M. H. Devoret, Phys. Rev. Lett. 76 , 3806 (1996).
- [36] B. T. Walker, J. D. Rodrigues, H. S. Dhar, R. F. Oulton, F. Mintert, and R. A. Nyman, arXiv:1908.05568.
Fahri Emre Ozturk Tim Lappe Goran Hellmann Julian Schmitt Jan Klaers Frank Vewinger Johann Kroha Martin Weitz
Fluctuation dynamics of an open photon Bose-Einstein condensate
Bosonic gases coupled to a particle reservoir have proven to support a regime
of operation where Bose-Einstein condensation coexists with unusually large
particle-number fluctuations. Experimentally, this situation has been realized
with two-dimensional photon gases in a dye-filled optical microcavity. Here, we
investigate theoretically and experimentally the open-system dynamics of a
grand canonical Bose-Einstein condensate of photons. We identify a regime with
temporal oscillations of the second-order coherence function $g^{(2)}(\tau)$,
even though the energy spectrum closely matches the predictions for an
equilibrium Bose-Einstein distribution and the system is operated deeply in the
regime of weak light-matter coupling. The observed temporal oscillations are
attributed to the nonlinear, weakly driven-dissipative nature of the system
which leads to time-reversal symmetry breaking.