Observation of nonlinear response and Onsager regression in a photon Bose-Einstein condensate

The quantum regression theorem states that the correlations of a system at two different times are governed by the same equations of motion as the single-time averages. This provides a powerful framework for the investigation of the intrinsic microscopic behaviour of physical systems by studying their macroscopic response to a controlled external perturbation. Here we experimentally demonstrate that the two-time particle number correlations in a photon Bose-Einstein condensate inside a dye-filled microcavity exhibit the same dynamics as the response of the condensate to a sudden perturbation of the dye molecule bath. This confirms the regression theorem for a quantum gas, and, moreover, demonstrates it in an unconventional form where the perturbation acts on the bath and only the condensate response is monitored. For strong perturbations, we observe nonlinear relaxation dynamics which our microscopic theory relates to the equilibrium fluctuations, thereby extending the regression theorem beyond the regime of linear response.


INTRODUCTION
The application of linear response theory to systems that are subject to perturbations lies at the heart of many fundamental phenomena in physics [1], including electromagnetic wave propagation in optical media, structure factors in condensed matter systems, or superfluid phases in quantum gases [2][3][4].For strong perturbations, the extension of this concept to nonlinear response has facilitated our understanding of ubiquitous effects such as higher-harmonic generation in optics [5] or the emergence of turbulent flow in coldatom and exciton-polariton systems driven far from equilibrium [6,7].When considering equilibrium systems, the linear response behaviour is commonly governed by the fluctuationdissipation theorem [8], which states that the intrinsic fluctuations of a system are connected to the absorptive part of a response function by thermal energy.The Onsager-Lax theorem, covering situations where the magnitude of the fluctuations is small, remains valid also for systems out of equilibrium and describes a universal relationship between the correlations and the response dynamics, as has been theoretically shown for irreversible processes in classical and quantum systems [9][10][11][12][13][14][15][16][17].The usual regression theorem, which states that the two-time averages ⟨A(t)B(0)⟩ of two observables A and B obey the same equations as the one-time averages ⟨A(t)⟩ for Markovian systems [12], links the system's fluctuations to the linear response, but more recent works have theoretically addressed nonlinear response as well [18].Experimentally, it has been indirectly verified by measurements of Onsager's reciprocity relations in the classical domain, e.g., in thermoelectric systems or thin films [19,20].A direct experimental test of this central theorem of statistical physics, however, by independent measurements of the temporal response to a perturbation and of the system's fluctuation dynamics has so far not been carried out for quantum gas systems [21][22][23].
To examine the regression theorem for a quantum gas, we investigate the reservoir-induced dynamics of a Bose-Einstein condensate (BEC) of photons in a dye-filled optical microcavity.In the used experimental platform, a twodimensional photon gas is coupled radiatively to a reservoir of dye molecules [24][25][26].Previous work has, by observing the statistical number fluctuations that result from an effective particle exchange with the molecular reservoir, reported a non-Hermitian phase transition [27,28], verified a fluctuationdissipation relation connected to a reactive response function [29], and realised protocols to temporally perturb the reservoir [30,31].Moreover, theory work on this system has employed the regression theorem to calculate two-point correlations from the relaxation dynamics [32][33][34].
Figure 1a illustrates how a photon BEC couples to a reservoir of dye molecules, providing a benchmarking platform for the regression theorem in complex many-body quantum systems of both matter and light.In a unique fashion, the system allows one to connect the linear response dynamics after a weak perturbation (Fig. 1a, left panel) with the intrinsic fluctuations driven by coloured noise (Fig. 1a, middle) and with the nonlinear response after a strong perturbation (Fig. 1a, right).Unlike the more usually encountered situation for the quantum regression theorem where a system (here the BEC) itself is perturbed, the perturbation here acts on the reservoir part (molecules), the excitation is transferred to the system by light-matter interactions (illustrated by a spring), and only then the photon number dynamics is witnessed by an observer before it is damped out by dissipation of photons to the environment.Despite its microscopic complexity, at the mean-field level the photon-molecule system can effectively be described as a damped harmonic oscillator, which exhibits a pronounced nonlinearity for large perturbations.
In this study, we measure the nonequilibrium response dynamics of a photon Bose-Einstein condensate after a sudden perturbation of its equilibrium dye reservoir.By comparing the response to the independently measured number fluctua- coupled to a molecule reservoir (bottom, red), forming a quantum dissipative oscillator.A weak sudden perturbation of the reservoir (left) leads to a linear response of the average single-time photon number ⟨n(t)⟩.Immersion into a "coloured-noise" bath (middle) leads to number fluctuations δn, or alternatively, two-time correlations ⟨δn(t)δn(0)⟩, for which the regression theorem predicts the same dynamics.For a strong perturbation (right), the response becomes nonlinear.b, Dye-filled optical microcavity with photon Bose-Einstein condensate, which is perturbed by a laser pulse irradiated on the dye reservoir (inset).The fluctuation and response dynamics are recorded on a photomultiplier (PMT).c, In a steady state, when cavity losses are compensated by pumping the dye, photons (blue) and molecules in the ground and excited state (red) are in equilibrium due to absorption and emission events.After a sudden perturbation of the excited molecules, also the condensate population is perturbed and both systems relax back to their stationary values.d, Left: Temporal evolution of the condensate averaged over several time traces.Right: Zoom-in on time window used to analyse the photon number response dynamics (top row), and time window used for the fluctuation analysis (bottom).Left panels show single-shot data and right panels the averaged data.
tions of the condensate, we first experimentally confirm the validity of the regression theorem for the optical quantum gas in the limit of weak perturbations.Further, for stronger perturbations we observe the emergence of nonlinear dynamics that is attributed to the saturation of the particle reservoir.Notably, the seemingly violated regression theorem is restored by a theoretical model that captures the nonlinear response dynamics, whose relevant parameters are the same as in the linear response.Such a nonlinearity in the condensatebath system, which has been theoretically predicted also for exciton-polariton systems [35], forms the basis for future studies into the properties of elementary and topological excitations within lattices of photon condensates [36][37][38].

Experimental scheme
Our photon Bose-Einstein condensates are prepared inside an optical microcavity filled with a dye solution of refractive index ñ ≈ 1.44 and concentration 1 mmol L −1 , see Fig. 1b; for details see Ref. [32] and Methods.The microcavity is formed by two spherical mirrors with a reflectivity > 99.998% and 1 m radius of curvature.At the used microcavity length D 0 ≈ 1.5 µm, the free spectral range of the cavity becomes as large as the emission and absorption spectral profiles of the dye medium, which restricts the photon dynamics to the two transverse degrees of freedom at a fixed longitudinal mode number q = 7. Inside the microcavity, the photons behave as a two-dimensional gas of bosons with effective mass m ph = πℏqñ/(D 0 c) ≈ 10 −35 kg and quadratic dispersion; the minimum photon energy m ph (c/ñ) 2 = ℏω c ≈ 2.1 eV is given by the energy of the transverse ground mode.In addition, the mirror curvature induces a harmonic trapping potential of frequency Ω/(2π) ≈ 40 GHz for the photons.The emission and absorption rates B em (ω) and B abs (ω) of the dye medium at T = 300 K fulfil the Kennard-Stepanov relation B em /B abs ∝ exp(−ℏ∆/k B T ), which depends on the detuning ∆ = ω − ω zpl of the photon frequency from the zero-phonon line ω zpl [27].By absorption-emission cycles with dye molecules (10 −12 s timescale), the photons thermalise to the rovibronic temperature of the dye before leaving the cavity (10 −9 s).Above the critical photon number 2 ≈ 80000, the photon gas exhibits Bose-Einstein condensation [24][25][26].Despite the thermalisation mechanism, the photon Bose-Einstein condensate realises a weakly-dissipative macroscopic quantum system due to losses, e.g., from photons leaking through the cavity mirrors.To establish a steady state of photons and molecules, the dye medium is externally pumped.

Theoretical description
Under steady-state conditions, the weakly-dissipative character of the open condensate system is evident from an exceptional point in the second-order temporal correlations [28,34].The underlying stochastic number fluctuations are caused by an effective particle exchange between the condensate and a large reservoir of excited molecules [27].In the present work, we access the open-system dynamics of the condensate by investigating the photon number response after a controlled sudden perturbation of this reservoir.This allows us to test the validity of the regression theorem for the optical quantum gas, as well as to explore it beyond the regime of linear response.To begin with, we derive an analytical expression for the response dynamics [39] from two coupled rate equations for the number of photons n and excitations X = n + M ↑ , respectively, given by dn/dt = B em M ↑ (1+n)−B abs M ↓ n−κn and dX/dt = P M ↓ −κn−Γ sp M ↑ (Methods).Here M ↓,↑ denotes the number of molecules in their ground (↓) and electronically excited (↑) states, P the pump rate, κ the cavity loss rate, and Γ sp the spontaneous decay rate to unconfined modes.A steady state is established if the losses are compensated for by a constant pumping of rate P = κ⟨n⟩/⟨M ↓ ⟩ + Γ sp ⟨M ↑ ⟩/⟨M ↓ ⟩, where ⟨...⟩ denotes temporal averaging, see Fig. 1c.
A time-dependent perturbation P (t) of the molecule reservoir, however, drives the photon condensate away from the steady state.The resulting photon number evolution n(t) = ⟨n⟩ exp{B em [1+exp(ℏ∆/k B T )] t 0 m(t ′ )dt ′ } depends on the strength of the perturbation and therefore on the deviation m(t) of the excited molecule number from the steady state number that would correspond to the instantaneous photon number n(t) (see Methods).The implicit interplay between m(t) and n(t) also implies that a subsequent variation of the photon number leads to a change of the number of excited molecules.As detailed in the methods, approximating the time integral over the perturbed molecules yields a nonlinear expression for the photon number response R(t) = n(t)−⟨n⟩: Here m 0 = m(0) gives the number of molecules excited by the pulse perturbation at t = 0.The quantities s ± = −δ ± δ 2 − ω 2 0 depend on the system parameters ∆, B em , T, Γ sp , κ, molecule number M , and the steady-state photon number ⟨n⟩, where the damping rate δ and the oscillation frequency ω 0 determine whether a biexponentially damped (δ > ω 0 ) or an oscillatory (δ < ω 0 ) response is expected.Note that s ± are the eigenvalues of the non-Hermitian matrix which describes the linearised equations of motion for the variation of the photon number n(t) − ⟨n⟩ and excitation number X(t) − ⟨X⟩; for a discussion see the Methods and Ref. [28].Expanding eq. ( 1) yields a linear expression for the response R(t) ≃ ⟨n⟩B em [1 + exp(ℏ∆/k B T )]m 0 (e s+t − e s−t )/(s + − s − ).Both expressions are used to analyse the response dynamics, but eq. ( 1) describes the measured response more accurately in the regime of strong perturbations, while the linearisation is valid in the limit of small m 0 .The expected second-order coherence at time delay τ on the other hand can, by utilizing the regression theorem, be written as g (2) (τ ) ∝ (s + e s+τ − s − e s−τ )/(s + − s − ) ∝ dR(τ )/dτ (see Methods).Physically, the derivative follows from the fact that it takes some time for the molecules excited by the laser pulse to be converted into photons.This gives a delay for the response function with respect to the correlation function, in contrast to a direct perturbation of the photon number.Unlike the response, the fluctuation dynamics g (2) (τ ) is intrinsically linear due to the grand canonical condition for large fluctuations being realised only for large reservoirs which basically cannot be saturated.

Experimental protocol
To probe the response and fluctuation dynamics, we first prepare a steady-state photon BEC by quasi-cw optical pumping of the dye molecule reservoir over 600 ns at 532 nm wavelength, see Figs. 1b,c.After roughly 200 ns, a short laser pulse of 28 ps duration (also at 532 nm) irradiates the microcavity to perturb the reservoir.Part of the emission leaking from the microcavity is filtered in momentum space and for polarisation, and the photon number evolution in the transmitted condensate mode is recorded using a photomultiplier (Methods).Figure 1d shows an example of the condensate evolution after averaging over many time traces, from which the response dynamics is visible in a time window of 20 ns after the perturbation.The fluctuation dynamics is determined by analysing the second-order correlations g (2) (τ ) from individual time traces.Throughout all measurements the system parameters B em = 25 kHz, ℏ∆ = −3.87kB T, Γ sp = 200 MHz remain fixed, while the total molecule number M = 5.4(15) • 10 9 , and mirror transmission κ = 6.4 (10) GHz are obtained from fits to the data (Methods).

Second-order correlation and response dynamics
Figure 2 shows measured second-order correlation functions g (2) (τ ) and photon number responses R(t) along with fits for two steady-state photon numbers ⟨n⟩, realised by varying the quasi-cw pump power.We are first interested in the linear response of the photon condensate, so we use only relatively small perturbation pulse powers that weakly change Frequency, biexp. the condensate population by δn(t max )/⟨n⟩ = 0.31 (7) on average, where t max denotes the time when the photon population has reached its maximum.Accordingly, the response data is fitted using the linearised form of R(t) discussed above.While generically g (2) (τ ) = 1 and R(t) = 0 are found for large τ and t, both data sets display distinct dynamics.For small ⟨n⟩ = 1740, both g (2) (τ ) and R(t) decay biexponentially, while for large ⟨n⟩ = 14250 a damped oscillation of the fluctuations and the response is observed.The quantitative agreement of the dynamics for each photon number is seen in the respective values of the parameters δ, ω 0 , see Fig. 2c, which determine the eigenvalues s ± .This gives evidence that the intrinsic number fluctuations and the response dynamics of the photon condensate to an external perturba-tion of the reservoir are governed by the same microscopic physics.Qualitatively, the agreement is clearly visible when forming the derivative of the fitted response dR(t)/dt, see the inset of Fig. 2b, which well resembles the corresponding fit for g (2) (τ ) from Fig. 2a.Note that the biexponential and oscillatory dynamics are distinctly characterised by real and complex-valued s ± , respectively, which allows to identify a transition point between both dynamical regimes at the degeneracy s + = s − , known as an exceptional point [28].

Regression theorem in linear regime for weak perturbations
To systematically verify the universal relationship between R(t) and g (2) (τ ) in the linear response regime, we next study the eigenvalues s ± of the condensate dynamics as a function of the steady-state population ⟨n⟩. Figure 3a shows the measured damping rate Re(s ± ) and oscillation frequency Im(s ± )/(2π) for the response and fluctuation dynamics (symbols), respectively, along with the eigenvalue prediction from eq. ( 1) (lines).Note that for the response R(t) we only show data for the smallest experimentally realised perturbation powers (as in Fig. 2), where the linear model is expected to be well applicable.At photon numbers below the one at the exceptional point ⟨n⟩ EP ≈ 2000, two branches of damping constants with vanishing oscillation frequency indicate the biexponential regime, while for larger photon numbers the observed merging of Re(s ± ) and bifurcation in Im(s ± ) highlight the regime of oscillatory dynamics.
Figure 3a can be understood as a nonequilibrium phase diagram of the photon dynamics, where ⟨n⟩ presents a control parameter to tune between both phases.The associated opening of a gap in the complex plane along the imaginary axis at the phase transition is visible in the complex eigenvalue spectra of g (2) (τ ) and R(t) in Fig. 3b; symbol colours indicate the photon number, which parametrises the spectral trajectory of s ± (⟨n⟩).As the photon number is increased, the real-valued s ± (biexponential dynamics) move from {Re(s + ), Im(s + )} = {−∞, 0} and {Re(s − ), Im(s − )} = {0, 0} along the real axis and coalesce at the exceptional point {Re(s ± ), Im(s ± )} ≈ {−0.8, 0}ns −1 .As ⟨n⟩ is increased further, the eigenvalues s ± separate again and move into the imaginary plane (oscillatory dynamics).
The agreement between the measured linear response, the fluctuation dynamics, and the theory prediction provides an experimental benchmark for the regression theorem in optical quantum gases.Direct evidence for this conclusion is given in Fig. 3c, which compares the eigenvalues s ± of the response measurement to the values of s ± obtained from the fluctuation measurement for pairwise corresponding photon numbers.In order to represent the two sets (s + and s − ) of complex-valued data, we show the four contributions Re(s + ), Re(s − ), Im(s + ) and Im(s − ).We find excellent agreement between both sets of eigenvalues, as evident from the data aligning with a linear curve of slope one (black line).Interestingly, for large condensate populations currently not accessible in the experi- ment and shown in Fig. 3d, our theoretical model for the photon dynamics predicts a second exceptional point where again s + = s − , or equivalently δ = ω 0 (see Methods for the corresponding functions).Physically, the re-emerging biexponential phase results from the different asymptotic scalings of δ ∼ ⟨n⟩ and ω 0 ∼ ⟨n⟩, such that oscillations are damped out not only at average photon numbers smaller that ⟨n⟩ EP but also in the limit of very large ⟨n⟩.

Nonlinear response for strong perturbations
We next generalise our study of the regression theorem to the nonlinear regime by successively increasing the reservoir perturbation strength m 0 induced by the pulse laser.Figure 4a shows complex eigenvalue spectra s ± (⟨n⟩) obtained from fitting either the linear (top row) or nonlinear (bottom) expression of R(t) to the recorded condensate time traces.While for the weak perturbation the linear and nonlinear analysis yield similar results, a deviation from linear response theory is observed for the larger perturbations, as highlighted by the orange shading.In contrast, fitting the nonlinear expression in eq. ( 1) restores the eigenvalue spectrum of the response dynamics to agree well with the theoretical prediction for the fluctuation dynamics.This improvement demonstrates that the strongly perturbed photon condensate exhibits a pronounced nonlinearity, which occurs in our system when the molecule reservoir is saturated: almost all the excitations that are introduced by the perturbation pulse are then converted into photons leading to a large relative variation of the photon number, that activates the nonlinearity in the stimulated emission and absorption dynamics.Figure 4b shows the residuals i (R exp,i − R fit,i ) 2 /N /⟨n⟩ of the linear and nonlinear fit as the initial pulse perturbation is increased.Here, N denotes the number of recorded samples.Both the linear and the nonlinear fit yield similar and small residuals for weak perturbations below m 0 ≈ 5 • 10 4 (extracted from the fit), confirming that the condensate-reservoir coupling can here be well described by linear dynamics.Beyond this value, however, the residuals exhibit a splitting.While the linear fit performs significantly worse (i.e., the residuals grow), the nonlinear fit residuals remain close to their initial values, demonstrating that the nonlinear expression is more accurate in describing the response dynamics of the photon condensate in the case of relatively strong perturbations of the molecule reservoir.Despite the improved accuracy of the nonlinear model, we note that for even stronger perturbations, as theoretically expected, also our nonlinear description gradually loses its accuracy due to a linearisation in m (Methods).to the molecule reservoir is directly visible when we compare the measured oscillatory response R(t) in Fig. 4c for a weak and a strong perturbation strength with m 0 ≈ 2.7 • 10 4 and m 0 ≈ 9.3 • 10 4 , respectively.The solid lines show fits based on the linearised and full nonlinear expressions.For the weak perturbation, we find both models to describe the experimental data equally well.For the larger perturbation strength, the measurement is more accurately fitted by the nonlinear expression, as highlighted in the zoom-in view on the 5 to 12 ns time range in the inset of Fig. 4c.

DISCUSSION
To conclude, we have measured the response dynamics of a photon Bose-Einstein condensate after a sudden perturbation of the reservoir, which before the perturbation forms an equi-librium steady-state with the condensate.Comparing the response dynamics with the number fluctuations has enabled the experimental verification of the regression theorem for optical quantum gases.Specifically, we have identified a nonlinear response of the BEC at strong driving.Here an extended regression theorem of the form g (2) (τ ) − 1 ∝ d/dt ln[1 + R(t)/⟨n⟩] holds, which in the limit of small perturbations agrees with the conventional linear form.For the future, the demonstrated scheme to perturb photon condensates constitutes a novel tool for exploring Kibble-Zurek dynamics [40,41], or vortex turbulence [7,36] in optical quantum gases in tailored potentials [42,43].An extension of the reported single BEC-bath oscillator system to arrays of coherently coupled condensates linked with local reservoirs will enable the exploration of reservoir-induced transport dynamics and open-system topological states [37,38].Our findings also pave the way for studies of the time-dependent fluctuation-dissipation relation in photon condensates, which so far have only been confirmed for static reactive response functions [29].

METHODS Experimental Methods and Calibrations
For the steady-state excitation of the photon Bose-Einstein condensate inside the dye-filled microcavity, we use a cw laser at 532 nm wavelength (Coherent Verdi 15G).To minimise bleaching of the dye medium (Rhodamine 6G solved in ethylene glycol, concentration 1 mmol L −1 and zero-phonon line ω zpl = 2π × 550 THz) as well as pumping to nonradiative triplet states, the pump beam emission is temporally chopped into pulses of 600 ns at a repetition rate of 50 Hz by acoustooptic modulators.During the steady state, a mode-locked laser pulse of 28 ps duration at 532 nm wavelength (EKSPLA PL 2201) is irradiated onto the dye-filled cavity at the same repetition rate, i.e., there is one perturbation pulse during the 600 ns-long steady-state.The laser pulse instantaneously perturbs the number of excited dye molecules, and consequently also drives the photon condensate out of its steady state.The residual condensate emission transmitted through both cavity mirrors is used for the analysis of the experiment.On one side of the cavity, a microscope objective (Mitutoyo MY10X-803) collects the emission to measure the spectral distribution and the spatial intensity distribution of the photon gas.The cavity length is actively stabilised by monitoring the cutoff wavelength behind an Echelle grating on an EMCCD camera (Andor iXon 897).On the opposite cavity side, the cavity emission is filtered by truncating high-momentum states of the photon gas using an iris in the far field.The filtered condensate mode is detected using a photomultiplier (PHOTEK PMT 210) with a temporal resolution of 150 ps FWHM that is sampled by an oscilloscope (Tektronix DPO7354C) operated at 20 GSa/s sampling rate with a 2 GHz bandwidth.To calibrate the condensate population ⟨n⟩ against the recorded PMT voltage, the photon gas spectrum is fitted with a Bose-Einstein distribution [24].
To exclude systematic sources of errors in our PMTbased measurements of the temporal second-order correlations g (2) (τ ), we perform a benchmark with a HeNe laser at 632.8 nm wavelength, see Supplementary Fig. 1.For the coherent source, one expects g (2) (τ ) = 1 for all time delays.However, the obtained signal shows g (2) (τ ) > 1 up to τ ≈ 3 ns, an observation which we attribute to electronic noise in the detection system.To avoid a misinterpretation of the data arising from this artefact, we exclude data points at τ ≤ 3 ns from the analysis of the correlation dynamics.Moreover, radiofrequency noise collected by our detection system (e.g., during pulse picking in the pulse laser system) is visible in the averaged time traces of the photon condensate response at the time of the pulse emission.To mitigate this, an optical delay line has been implemented to shift the pulse arrival time at the microcavity with respect to the noise signal.The delay line is realised by a cavity of 24 m length (corresponding to a time delay of 80 ns), which is traversed by the pulse 6 times.To that end, all response measurements can be performed in a temporal region free of residual radiofrequencyinduced noise.

Theoretical Model
Using the Kennard-Stepanov relation for the absorption and emission rates of the dye medium, the rate equation for the number of photons n can be rewritten as We represent the of excited molecules as where is the number of excited molecules, which corresponds to the steady state with the average photon number equal to n.Then for the experimentally relevant case n ≫ 1, eq. ( 2) takes the form which leads to the following nonlinear relation between m(t) and the corresponding evolution of the photon number, starting from its value ⟨n⟩ at t = 0: To describe the dynamics of m(t) initiated by a sudden perturbation of the excited molecule number at t = 0, we utilise the rate equation for where ⟨M ↑ ⟩ ≡ M ↑,⟨n⟩ and with Inserting eqns.( 4), ( 6) and ( 8) into eq.( 7), one obtains a rather cumbersome nonlinear differential equation for t 0 m dt ′ , which, in general, cannot be solved analytically.To proceed further, we assume that m is relatively small and hence it can be estimated from the linearised version of the aforementioned equation with where is the photon number relaxation rate [32].Equation (10) has the solution with Here, m 0 = m(0) is the number of molecules excited by the initial perturbation.Inserting eq. ( 14) into (6), we can express the time dependence of the photon number in an explicit form: Relation between R(t) and g (2) (t) The regression theorem implies that the response dynamics R(t) = n(t) − ⟨n⟩ of the photon condensate after a perturbation of the molecule reservoir is related to the condensate's second-order coherence g (2) (t) = ⟨n(t)n(0)⟩/⟨n⟩ 2 by g (2) (t) − 1 ∝ dR(t)/dt.In linear response, for small deviations around the mean photon number n = ⟨n⟩ + ∆n and excitation number X = ⟨X⟩ + ∆X, we have the following equations of motion for the time evolution after the system has been perturbed [28]: with the steady state photon variance ⟨δn 2 ⟩ = M eff ⟨n⟩ 2 /(M eff + ⟨n⟩ 2 ).The eigenvalues of this matrix correspond to s ± from eq. (15).Note that we here extend the description of ref. [28] by including losses from the spontaneous decay of the molecules into unconfined modes, which improves the theory description of the experimental observation.The first equation of the matrix form in eq. ( 17) expresses (i) that deviations in the photon density at constant X relax at the rate Γ and (ii) that changes in X lead to a change in the photon number.The second equation of the matrix form expresses that excitations are lost through cavity losses and spontaneous emission.
Here δn(t|δn 0 ) is the average photon deviation starting from a deviation δn 0 at time t = 0. From eqns.(17) and for ∆n = δn 0 , ∆X = 0 one finds From eqns.( 19) and ( 20), we then obtain where we have set δn 2 0 to the time averaged photon variance ⟨δn 2 ⟩.This expression contains the same rates (and frequencies) as the response function in eq. ( 18) but does not show identical time dependence.From the derivation of the correlation function, one can see that a perturbation in the number of photons at constant X would give a density response with the same time dependence as G (2) (t).Such a perturbation is hard to implement experimentally.By comparing eqns.(18) and ( 21) one finds Using G (2) (t) = [g (2) (t) − 1]⟨n⟩ 2 , this corresponds to the relation stated in the beginning of this section.In the case of an oscillating density response, the derivative implies that the correlation and response functions have a phase shift of π/2.Physically, this can be understood from the fact that an external laser pulse excites the molecules, which take some time to be converted into photons.This gives a delay for the response function with respect to the correlation function.In Fig. 2 one sees that the correlations show a minimum as a function of time while the response does not.The correlation function G (2) (t) must become negative as can be seen from the differential relation in eq. ( 22) and ∆n(t → ∞) = 0 in the presence of losses.Physically, it is a consequence of a positive photon number fluctuation at some moment to imply larger losses, which reduces the expected photon number later on.

Linearity of Fluctuation Dynamics
The particle number fluctuations and the corresponding second-order correlations of a photon Bose-Einstein condensate are governed by linear dynamics even under grand canonical statistical conditions.Here a simple analytical reasoning for this statement is presented.The rate equations for the coupled photon-dye system can be employed to identify two competing rates, which determine whether nonlinear effects in the photon dynamics can or must not be neglected [32].In the lossless case (κ = Γ sp = 0), the rate equation for a photon number fluctuation away from its steady-state value reads where Γ (⟨n⟩) is defined in (13) (we indicate here explicitly its photon number dependence).The effective reservoir size M 0 eff is given by eq. ( 9) with κ = 0. Nonlinear effects become relevant for B em [1 + exp(ℏ∆/k B T )] ∆n ∼ Γ (⟨n⟩).Inserting δn (see Section 'Relation of R(t) and g (2) (t)'), one has Supplementary Fig. 2 shows both rates as a function of ⟨n⟩ for different dye-cavity detunings, which changes the effective reservoir size.For all curves the second-order correlation rate (r.h.s) exceeds the nonlinear term (l.h.s.), meaning that a photon BEC driven by reservoir-induced fluctuations to good approximation always exhibits linear dynamics.Our experimental data in Figs. 2, 3 confirm this prediction, showing that the second-order correlation dynamics are well described by the derivative of the linear expansion of eq. ( 1).

Fitting of Experimental Data
To analyse the dynamics of the two-time correlations, we fit the second-order correlation data with g (2) (τ ) = a s + e s+(τ +∆τ ) + s − e s−(τ +∆τ ) s + − s − + b, (25) where the eigenvalues of the dynamics are given by s ± = −δ ± δ 2 − ω 2 0 .Accordingly, for the response function in the nonlinear form we use the fit function R nl (t) = a exp b e s+(t+∆τ ) + e s−(t+∆τ ) s + − s − − a (26) and for the linear form we fit R lin (t) = a ′ e s+(t+∆τ ) + e s−(t+∆τ ) s + − s − The fit parameters δ and ω 0 resemble the damping constant and natural frequency of a harmonic oscillator; the time delay∆τ together with a, a ′ , b are treated as free parameters for each fit.From the analysis of the condensate response dynamics, we find that both the linear and nonlinear fit can be used to describe the experimental data depending on the strength of the perturbation pulse, as shown in Figs.4b,c of the main text and in Supplementary Fig. 3.For a quantitative comparison of the measured eigenvalues s ± with theory, several system parameters are required: On the one hand, the spontaneous molecule decay rate to unconfined optical modes Γ sp ≈ 200 MHz and the Einstein coefficient for emission B em = 25 kHz at the dye-cavity detuning ℏ∆/k B T = −3.87(corresponding to a cutoff wavelength λ c = 570 nm) are known [28,32].The molecule number M and the cavity loss rate κ, on the other hand, need to be determined.For this, we compare our experimental results to a full numerical solution of the coupled photon-molecule rate equations and minimise the difference by varying the two parameters.To achieve this in a consistent way, all recorded time traces (i.e., for all perturbation powers and all average photon numbers) are fitted simultaneously.This numerical approach is motivated by the fact that even the nonlinear expression for the response dynamics in eq. ( 1) is not exact, because in its derivation m was assumed to be small.We obtain a molecule number M = 5.4(15) • 10 9 and a cavity loss rate κ = 6.4(10)GHz.The results are consistent with the corresponding values from analysing the secondorder correlation function g (2) (τ ), M = 6.6(7) • 10 9 and κ = 6.6(10)GHz, which are obtained from fits of the linear expression in eq. ( 25).As discussed above, this is justified because the fluctuation dynamics are expected to obey linear equations.
Finally, the numerical data allows us to cross-validate the fit results of the linear and nonlinear expressions to the experimental data.Supplementary Fig. 3b shows a map of linear and nonlinear residuals obtained from fitting numerically calculated data as a function of the photon number ⟨n⟩ and perturbation strength m 0 .Fitting the nonlinear expression always yields smaller residuals compared to the linear fit at the respective coordinates; Supplementary Fig. 3c shows the residuals versus m 0 when averaged for average photon numbers ⟨n⟩ ≥ 4000, which qualitatively agrees with the experimental observation of Fig. 4b that the nonlinear model can describe the data more accurately at large m 0 .

Fig. 1 .
Fig.1.Experimental scheme to probe the regression theorem.a, Mechanical analogue of a photon Bose-Einstein condensate (top, blue) coupled to a molecule reservoir (bottom, red), forming a quantum dissipative oscillator.A weak sudden perturbation of the reservoir (left) leads to a linear response of the average single-time photon number ⟨n(t)⟩.Immersion into a "coloured-noise" bath (middle) leads to number fluctuations δn, or alternatively, two-time correlations ⟨δn(t)δn(0)⟩, for which the regression theorem predicts the same dynamics.For a strong perturbation (right), the response becomes nonlinear.b, Dye-filled optical microcavity with photon Bose-Einstein condensate, which is perturbed by a laser pulse irradiated on the dye reservoir (inset).The fluctuation and response dynamics are recorded on a photomultiplier (PMT).c, In a steady state, when cavity losses are compensated by pumping the dye, photons (blue) and molecules in the ground and excited state (red) are in equilibrium due to absorption and emission events.After a sudden perturbation of the excited molecules, also the condensate population is perturbed and both systems relax back to their stationary values.d, Left: Temporal evolution of the condensate averaged over several time traces.Right: Zoom-in on time window used to analyse the photon number response dynamics (top row), and time window used for the fluctuation analysis (bottom).Left panels show single-shot data and right panels the averaged data.

Fig. 4 .
Fig. 4. Nonlinear BEC response.a, Complex eigenvalue spectra for increased reservoir perturbations with δn(tmax)/⟨n⟩ = {0.59(10),1.15(23), 1.25(25)} and m0 = {5(2), 9(1), 11(3)} • 10 4 .Top row: Eigenvalues s± (data points) obtained from a linear fit gradually deviate from theory (solid lines) due to a saturationinduced nonlinearity (shaded area).Bottom: s± extracted from the nonlinear fit restore the agreement with theory, demonstrating a beyond-linear regression relation.Diamond symbols correspond to the s+ branch (red line), circles to the s− branch (blue line).b, Residuals of linear and nonlinear fits versus increasing m0 averaged for ⟨n⟩ ≥ 4000 indicate the improved accuracy of the nonlinear fit.Dashed line gives noise floor.c, Photon response after a weak (left) and strong perturbation (right) for ⟨n⟩ ≈ 8250, along with linear (blue line) and nonlinear (red) fits.Error bars show standard fitting errors in a, and standard statistical errors in b.