Spatiotemporal coherence of non-equilibrium multimode photon condensates

We report on the observation of quantum coherence of Bose-Einstein condensed photons in an optically-pumped, dye-filled microcavity. We find that coherence is long-range in space and time above condensation threshold, but short-range below threshold, compatible with thermal-equilibrium theory. Far above threshold, the condensate is no longer at thermal equilibrium and is fragmented over non-degenerate, spatially overlapping modes. A microscopic theory including cavity loss, molecular structure and relaxation shows that this multimode condensation is similar to multimode lasing induced by imperfect gain clamping.

We report on the observation of quantum coherence of Bose-Einstein condensed photons in an optically-pumped, dye-filled microcavity. We find that coherence is long-range in space and time above condensation threshold, but short-range below threshold, compatible with thermalequilibrium theory. Far above threshold, the condensate is no longer at thermal equilibrium and is fragmented over non-degenerate, spatially overlapping modes. A microscopic theory including cavity loss, molecular structure and relaxation shows that this multimode condensation is similar to multimode lasing induced by imperfect gain clamping. Quantum condensation and coherence are intimately linked for ensembles of identical particles. Condensation, defined by a macroscopically large fraction of all particles being in a single state (usually the ground state [1,2]) is typically associated with coherence as seen in the firstorder correlation function, which is proportional to the visibility of fringes of an interference measurement [3].
While observation of thermal equilibrium and macroscopic occupancy of the ground state are sometimes considered sufficient proof of Bose-Einstein condensation (BEC), the enhancement of coherence brought by BEC means that interferometry is one of the most urgent measurements to be made with a condensate [4,5]. Where thermal equilibrium is not completely reached, coherence is the defining characteristic of non-BEC quantum condensation, e.g for semiconductor exciton-polaritons [6][7][8][9] and organic polaritons [10,11]. In non-ideal Bose gases, such as ultracold atoms, interactions tend to reduce but not destroy the coherence [12][13][14].
Photon condensates in dye-filled microcavities are weakly-interacting [15][16][17][18], inhomogeneous [19,20], dissipative Bose gases close to thermal equilibrium at room temperature [21][22][23][24][25]. It is worth noting that the physical system has some similarities to a dye laser, with the decisive difference being that lasing is necessarily a non-equilibrium effect whereas photons can also undergo BEC in thermal equilibrium. Consequently BEC implies macroscopic occupation of the ground state independently of the pump properties, whereas a laser is characterized by a large occupation of exactly the mode that is most strongly pumped [26].
Unique among physical realisations of BEC, in dyemicrocavity photon BEC the particles thermalise only with a bath and not directly among themselves. This implies that the establishment of phase coherence in the condensation process is necessarily mediated via indirect interactions through by the dye, i.e. a system whose fast relaxation renders all mediated interaction incoherent. Notably, coherence in photon condensates has not yet been systematically measured.
Condensates with macroscopic occupation of two or more states without phase relation are called fragmented [27]. Whereas strong, attractive interactions favour fragmentation, repulsive interactions stabilise a single condensate mode [2,28,29]. Nevertheless, fragmentation has been observed using ultracold atoms in multiple spin states [30], or separated spatial modes [31,32]. Fragmented, dissipative condensates with spatially separated states have been seen in polaritons in semiconductors [33,34], and organic solids [35]. It has been proposed that for driven, dissipative bosonic systems, multimode condensation is a general non-equilibrium phenomenon [36], when driving happens faster than dissipation (such as loss, thermal equilibration or spatial redistribution).
Below threshold pump power, P th , the coherence time T and length L of the thermalised light are expected to be of order h/k B T 0 0.15 ps and λ dB = hλ 0 c / 2πk B T 0 n 2 L 1.5 µm where λ 0 590 nm is the wavelength of the lowest-energy cavity mode, T 0 = 300 K the temperature, c the speed of light in free space and n L the refractive index of the solvent filling the cavity [37][38][39]. The coherence time is predicted to be much greater above threshold than below [25,39], increasing further as the number of particles in the condensate increases, and the coherence length is expected to be at least as large as the whole condensate [40]. Multimode condensation may occur and its effect on coherence is not predicted [20].
In this manuscript we present measurements of the coherence properties of thermalised photons with both time delays and position shifts between the two arms of an interferometer. We describe how the controls and outputs of our imaging interferometer correspond to the underlying first-order correlation function, g (1) (r, r , τ ), as a function of positions r and r and time delay τ . We characterise the coherence time and length of the photon condensate as a function of pump power. Below and just above threshold, the measurements are compatible with thermal equilibrium theory. Far above threshold, the condensate fragments into multiple, spatially-distinct but overlapping, non-degenerate modes accompanied by a decrease of both spatial and temporal coherence. We interpret multimode condensation as a non-equilibrium phenomenon similar to gain saturation in lasers.
Our experiment starts by pumping a fluorescent dye in a high-finesse microcavity [19,22] in quasi-continuous conditions. The pump spot was elliptical with a minor axis of typically 50-60 µm diameter, and we use the 8 th longitudinal mode of the cavity with a cutoff wavelength of 590 nm. These parameters are known to produce near thermal-equilibrium conditions [20]. The cavity photoluminescence is imaged to infinity, then split. Half is split again and imaged onto an auxiliary camera and a spectrometer whose spectral resolution, about 0.2 nm, is insufficient to resolve the bare cavity modes which are separated by 0.05 nm. The other half is sent to an imaging Mach-Zehnder interferometer, as shown in Fig. 1. Each of the two arms of the interferometer has a delay line: one controlled by a piezo for the fine motion to scan over a fringe, the other controlled by a motor for coarse motion. The horizontal axis, x, of the last adjustable mirror in one arm is controlled by a motor, whose motion is converted to a shift in position of the image at the camera. Both outputs of the interferometer are sent onto a camera through a single imaging optic, imaged to two separate locations on the sensor. There is a linear-polarising filter in front of the camera, which increases the visibility of fringes.
The camera records a spatially resolved intensity distribution. If one arm of the interferometer is blocked, this corresponds to the intensity I(r) emitted from the cavity, i.e. the spatial profile of the condensate photoluminescence. Since pumping and detection in this experiment are quasi-continuous, all processes are stationary. Temporal resolution comes in terms of the path delay of the interferometer. The detected interferometer signal depends on r, r and τ , where r = (x, y) is the position on the camera, r = (x + δ x , y), with the displacement δ x introduced by one arm of the interferometer and τ is the temporal delay corresponding to the path-length difference between the two arms of the interferometer. The Michelson visibility of fringes V is directly related to the coherence g (1) : We scan the piezo-controlled delay, typically acquiring a set of 41 images, while maintaining all other parameters fixed. The principal result of our data analysis (explained in detail in Supplementary Material) is a four-dimensional set of visibility data, V (x, y, δ x , τ ). For any value of x and y we can extract a characteristic coherence time T (or length L), with a fit, usually to a Gaussian, in τ (or δ x respectively). In Fig. 1 (bottom, from left to right) we see an image of the photon condensate I(x, y), an image of the visibility of the same condensate V (x, y) δx=0,τ =0 and an image of its coherence time, T (x, y) δx=0 . Since the images are overlapped (δ x = 0), the visibility V is exactly equal to the coherence g (1) . Fig. 2 is generated by choosing a single pixel x 0 , y 0 and measuring the visibility V (τ, δ x ) x0,y0 as a function of long-range delay and image shift. The results are shown just above and far below threshold. Under inspection, the differences between V and g (1) were not noticeable, so we have presented V . Coherence time and length are inferred from a two-dimensional Gaussian fit to the data. Far below threshold, the length and time scales of coherence, 4.5 µm and 0.2 ps, are slightly longer than the thermal scales. This overestimation in both space and time is explained by finite spatial resolution of around 3 µm (see Supplementary Material). Above threshold the measured length, 14 µm, of the condensate is comparable to the size of the condensate itself, implying that the whole condensate shares one phase, as expected. The measured coherence time of 10 ps is also large, limited by condensate emission frequency fluctuations on timescales equal to the time between images, 200 ms. The condensate emission frequency variations are dominated by the variation of the cavity length at the limits of our locking scheme, which has a bandwidth of 20 Hz and resolution equivalent to about 0.05 nm in cavity cutoff wavelength [41]. Fig. 3 depicts the spectrum, image I(r) and visibility image V (r) for various pump powers above threshold. The condensate peak in the spectrum broadens with increasing pump power and breaks up into multiple peaks, i.e. the condensate splits into multiple non-degenerate modes. Fig. 3 shows three peaks, but we have seen up to five distinct peaks in some experimental runs where we used reduced pump spot sizes to lower the pump threshold.
FIG. 3. Normalised photoluminescence spectrum (left column), normalised image (middle column) and visibility image (right column) for various pump powers (rows, as labelled on the graph) above threshold P th . The spectrum broadens and splits into multiple modes, the condensate broadens in space and the visibility image fragments at higher powers. A small pump spot (30 ± 10 µm) was used to reduce threshold pump power.
Along with this non-degenerate multimode behaviour, the condensate broadens and the measured visibility develops a fragmented structure. With a poorer spectrometer resolution, the spatial broadening could have been taken as an indication of repulsive interactions [22]. Since, however, there is no blue-shift in the spectrum apart from variations of the cavity length [41], Fig. 3 gives clear evidence of the condensation of non-interacting photons in several modes, rather than quantum depletion reducing the condensate fraction.
We now ask: how does coherence vary across threshold and in the multimode regime? Are the non-degenerate modes coherent with each other, and is the multimode behaviour a sign of the breakdown of thermal equilibrium?
We define the spatial coherence length as the scale of a Gaussian fit to the visibility, as a function of shift between two images V (δ x ), and measure it for various pump powers. V (δ x ) and a cut through the photoluminescence intensity I(x) are shown in Fig. 4 (top), for two pump powers, one far below and one just above threshold. In Fig. 4 (middle), we compare experiments to a thermal equilibrium theory without dissipation (see Supplementary Material). The theory is based on a series expansion of the correlation function [42,43] which agrees with exact calculations [44]. There are no free parameters below threshold (solid lines), but the scaling of the horizontal axis is imprecise above threshold (shown as dashed lines), as the number of photons varies non-linearly with pump power [25].
Within the theory's range of validity (P P th ) there is quantitative agreement with the experiment. Far below threshold, the coherence length is much shorter than the characteristic size of photoluminescence, limited only by imaging resolution. With increasing power around threshold, coherence length grows, as the width of the emitted light decreases. At even higher powers, when the system enters the multimode regime, the intensity increases but the coherence length decreases to around 6 µm (approximately the harmonic oscillator length scale), indicating that the multiple modes are incoherent. The condensate is only partially coherent, in contradiction to the dissipative, thermal-equilibrium prediction of Ref. [40].
In Fig. 4 (bottom), we show the coherence time. Far below threshold, the shortest measured coherence time is limited by spatial resolution and marginal undersampling of the data [45]. Above threshold, an upper bound for coherence time is set by the vibrations of the cavity [41]. Barely into the multimode regime, coherence time decreases, suggesting no coherence between modes, in agreement with the spatial coherence data. Even though thermal-equilibrium theory (black lines, solid below threshold, dashed above) does not describe the coherence time as accurately as spatial coherence and intensity, it captures qualitatively the increase of temporal coherence as the threshold pump power is reached.
Condensate fragmentation cannot be explained by thermal-equilibrium processes. We therefore need to invoke a microscopic model that takes into account spatially-inhomogeneous pumping, molecular relaxation via the thermal bath of solvent vibrations, spontaneous emission and cavity loss [20]. In Fig. 5 we show the results of the model (see Supplementary Material for more details). For computation efficiency, the model is restricted to one dimension. With increasing pump power, condensation occurs first in the lowest mode, then subsequently in higher modes (left panel). When one mode reaches threshold, it locally clamps the excited state population of dye molecules, but sufficient gain remains at the edges that more modes can reach threshold. The multimode regime is reached for the lowest pump powers for a pump spot which is large enough to overlap with several spatial modes of the bare resonator (right panel). It is possible to extract approximate values for coherence length and time from the same microscopic model, and we find good qualitative agreement with Fig. 4 in all regimes, below, near and far-above threshold. In conclusion, we have observed first-order coherence of thermalised photons in a dye-filled microcavity, below and above condensation threshold. Spatiotemporal correlations are longer-range for the condensed than noncondensed state, and show increases in range even below threshold, in broad agreement with thermal-equilibrium theory. Above threshold, multiple modes are seen which is a signal of non-equilibrium, driven, dissipative processes [36]. There is no coherence between modes. In this case, a microscopic model shows that the fragmentation can be explained using concepts from laser physics, i.e. imperfect gain clamping. By generating inhomogeneous, nonlinear gain and loss processes it may be possible to create equivalent states in other trapped condensates such as polaritons [9,46] or atoms [47]. It would be intriguing to know if higher-order particle-particle correlations occur between modes even in the absence of phase coherence, and how superfluidity manifests itself in multimode condensates.
During review of this manuscript we became aware of related work on phase coherence of photon condensates [48]. We thank Jonathan Keeling and Henk Stoof for inspiring discussions, and acknowledge financial support from the UK EPSRC (via fellowship EP/J017027/1 and the Controlled Quantum Dynamics CDT EP/L016524/1), and the ERC (via ODYCQUENT grant).

S1. DATA ACQUISITION AND ANALYSIS
Our experimental apparatus is almost identical to [S1]. We pump a fluorescent dye in a high-finesse microcavity in quasi-continuous conditions, using 500 ns pulses, which are much longer than any thermalisation or cavity loss time scales in the system. The pulse repetition rate is varied so that the product of pump power and repetition rate is kept constant, up to 2000 mW laser output power, where the repetition rate is 500 Hz. Our maximum laser output power is 2200 mW, and about 50% of this light makes it into the cavity (through modulators and the cavity mirror). The repetition rate variation means that we have acceptable signal-to-noise over a very large range of pump powers. Images are integrated over, typically, 50-2000 ms. The pump spot was elliptical with an aspect ratio not far from one, and a minor axis between 50 and 60 µm diameter.

A. Acquisition
The interferometer signal is observed using a colour camera. All colour values are converted to monochrome by summing red and green channels. For 590 nm (our typical working wavelength) the sensor of our camera (PointGrey Grasshopper GS3-U3-23S6C-C) is roughly equally sensitive in both red and green channels. Both output ports of the interferometer are directed to the same sensor. Two areas of the sensor are assigned as in-phase (P ) and quadrature (Q). For each set of data, there is one image taken with one arm of the interferometer blocked, to allow us to align P and Q images.
To measure the visibility, we scan the voltage of the piezo controlling one of the delay lines while maintaining all other parameters fixed, so that the fine-scale delay time which we call τ f varies over about 3 periods of oscillation of the light (6 fs). The piezo voltage is typically set to 41 values covering about 3 complete fringes. At each τ f we take an image. Coarse-scale delays, τ , up to 300 ps are achieved with the motor on the other delay line, with * Present address: Analytical Science Division, National Physical Laboratory, Hampton Road, Teddington, Middlesex TW11 0LW, UK † Correspondence to r.nyman@imperial.ac.uk a precision of around 2 fs. The total delay is τ f + τ , but we treat τ f and τ as independent, since they vary over very different magnitudes, |τ f | |τ |. We overlap the images vertically (in the y axis) but with a motorised mirror mount induce a shift in the horizontal axis, δ x . A full set of data for measuring the spatio-temporal coherence consists of images in (x, y) for a variety of τ f , τ and δ x . Thus we build a five-dimensional data set I(x, y, δ x , τ, τ f ). The notation includes the calibration of imaging magnification, so all co-ordinates presented are scaled to the intra-cavity co-ordinates.
A larger dataset is taken by varying control parameters of the photon condensate itself. In this manuscript, we only vary pump power P p . Whenever the pump power is changed, the exposure and gain of the camera as well as the spectrometer are adjusted automatically to maximise dynamic range.
We make the reasonable assumption that the coherence varies only slowly compared to the oscillations of the light. The microscopic description of dye-microcavity photon condensates requires no processes faster than about 100 fs [S2-S6], which is much longer than the 2 fs oscillation period of the light, justifying this assumption.

B. Analysis
When analysing data, each pair of P and Q images is aligned and binned if needed. Images are aligned by minimising the average-sum-square of differences between P and Q image values with respect to the shift of coordinates, taking only pixels which are present in both images after shifting. Knowing this optimised shift, we can be sure that a given pixel in P corresponds to the same pixel in Q.
The intensities are I P (x, y, δ x , τ, τ f ) and I Q (x, y, δ x , τ, τ f ) for in-phase and quadrature respectively. In the analysis of data sets where the control parameters of the condensate (e.g. pump power P p ) vary, an array of values for each control parameter is constructed, e.g I P (x, y, δ x , τ, τ f , P p ). The full sixdimensional data consists of as many 25 000 images, taking up to 30 GB of memory. Often 4-by-4 pixel blocks are combined to reduce computational effort in analysis. The major challenge is to visualise this data. This data set can be analysed and visualised in a number of ways, as shown in the main text, Fig. 1.

Visibility estimator
Extraction of the visibility over the set of fine delays is the most important processing we do on the data: I P and Q (x, y, δ x , τ, τ f , P p ) → V (x, y, δ x , τ, P p ). Our estimator for the visibility is based on a Fourier method and is robust against amplitude noise and converts phase and frequency noise to a reduction in visibility, unlike the Michelson visibility criterion. Amplitude noise is intrinsic to the photon condensate [S7]. Frequency noise is largely due to cavity length fluctuations. The result is, for each P p , a four-dimensional visibility V (x, y, δ x , τ ), which has been reduced over τ f . Even averaging over much longer than cavity lifetime, there are some parameters of our experiment which are not well controlled, and so, close to threshold, there are large variations in photoluminescence intensity. Background noise (from readout of the sensor or background light) may also affect the interferometer outputs. The Michelson criterion of visibility (ratio of the differences between and the sums of maxima and minima of signal) is not robust against these kinds of noise. Phase noise is also present. To negate the effect of intensity fluctuations, we use arc-tangent of the ratio of the two quadratures: In principle, it is important to subtract background signal and noise, but in practice we notice no effect of so doing. The Michelson estimator for visibility is not robust because it does not make use of all the data available, only the maxima and minima. Instead we use of an estimator based on the Fourier transform. Using the autocorrelation of the arc-tangent data, we find the approximate frequency of the interference pattern. We then calculate the Fourier component of the arc-tangent signal at several frequencies more closely spaced than the Nyquist criterion, since the data typically only cover about 3 cycles. The Fourier component at angular frequency ω f is: for a sample of arctan data and fine delay times The maximum amplitude gives the amplitude of fringes. The fringe amplitude is then normalised by half the sum of the data to give the visibility V , as would be expected for a clear sinusoid with an offset which gave Michelson visibility V .
We have tested the estimator using a model of noise in the system. We model the in-phase and quadrature signals, I P and I Q as I P = A 0 +Ã cos ωτ f +φ +R P + C (S3) Variables with tildes are random variables. Phase noisẽ φ and readout noises for each channel,R P andR Q are drawn from Gaussian distributions. We define an underlying amplitude A 0 and visibility V 0 , which control the offset C = A 0 V 0 / (1 − V 0 ). The angular frequency in piezo-controlled delay time units is ω. The total amplitude is the sum of the underlying amplitude and amplitude noise,Ã. Amplitude noise in photon condensates is known to vary from normal to skewed with a long tail for large values, so we drawÃ from a scaled Poisson distribution.
We tested four types of estimators against this noise model: the Michelson criterion, a root-mean-square test, sinusoid-fitting and the Fourier-based method described above. The Michelson criterion is not robust against noise, and the root-mean-square method requires scanning over an exact integer number of fringe periods. Least-squares fitting a sine wave gave fits were not robust, often picking local not global optima. Finally, the Fourier-based method was found to be robust (especially against amplitude noise when using the arc-tangent data) and low-noise, although we find that it is slightly biased.
For extreme values of underlying visibility (near zero or unity), the estimator is biased towards 0.5. The extent of the bias depends on the amplitude of the noises. For experimentally realistic parameters, from an underlying visibility of 1.0 our estimator gives about 0.8. If the detector is weakly saturated, then the maximum inferred visibility may be further reduced. We do not correct for this in our presented data, as we are not fully certain that our noise model is complete, nor do we know the exact parameters that match our experiment. Also, the uncertainties in the inferred visibilities are about as large as the correction required to compensate for the bias.

S2. UNCERTAINTIES
The data presented in the main body of the manuscript are presented without error bars. There are two major sources of uncertainty in our experiments: fluctuations of threshold, and of visibility. The fluctuations are slightly faster than the typical time to acquire a data set, 1-5 hours.

A. Variability near threshold
In Fig. S1, we vary the pump power and measure the output light intensity averaged over a small region at the center of the image. The pump spot was smaller than use in the main text, to allow us to reach pump powers far above threshold here. Camera exposure and gain are set automatically for each power to avoid saturation. Pulse repetition rate was held constant at 500 Hz for this specific data set. Exposure is at least 2 ms, i.e. at least one 500 ns quasi-CW pump pulse is always detected. Below threshold, the signal to noise is less than unity. Above S3 threshold, the intensity grows linearly with pump power. A bi-linear fit reveals that the pump power at threshold is 260 mW. The measured threshold varies from one experimental run to another. However, many points below threshold show a condensate. In the range 180-250 mW, there may be a condensate or not. We do not know the cause of this variability. It is not the variation of cavity cutoff wavelength (although that does fluctuate on 10minute time scales). It may be related to polymerisation of the dye, which we know occasionally forms clumps requiring cleaning of the cavity. Although the data might appear to show hysteresis (memory effects), the pump power was varied in a random order. In the main manuscript, where power is varied, we reject data which are below threshold when they are expected to be above it, and vice versa. We accept only data that follow an approximately monotonic increase in output intensity as a function of pump power.

B. Visibility variability
Standard-deviation errors are of limited use near threshold as there is no reason to believe that visibility measurements for a given set of parameters are drawn from a normal distribution. They are as likely to be drawn from a bimodal distribution, corresponding to below-and above-threshold behaviour. We have tried to ascertain the limiting uncertainties in visibility away from threshold, by measuring a large sample of visibilities as a function of delay above threshold: see Fig. S2.
By oversampling, we can build sub-samples and evaluate their standard deviations. The highest sub-sampleaveraged visibility measured is about 0.7, although we know that our estimator is low-biased for such large visibilities. The largest shot-to-shot uncertainty in visibility is 0.15. For lower average visibility, the uncertainty in the visibility is lower. For example, we can measure non-zero visibility of 0.04 with a signal-to-noise of unity in a single measurement. Our noise model produces a FIG. S2. Limiting uncertainty in the visibility, significantly above threshold. Left: Visibility variation with delay for overlapped images as a specific pixel. Grey dots are individual data points. Points with error bars are averages over 13 points centred on the marker, with error bar being the standard deviation of that sample. Right: Standard deviation of visibility as a function of mean visibility. Maximum visibility inferred is 0.7. Our estimator is low-biased for the largest visibilities, and high-biased for the smallest visibility values. similar pattern, i.e. maximum inferred visibility 0.7 and standard deviation > 0.1, only with unrealistically large phase noise, of amplitude at least π/4 radians. We conclude that there is intrinsic noise in the visibility which does not come from our measurement apparatus or visibility estimator.

S3. THERMAL EQUILIBRIUM THEORY OF BOSE GAS COHERENCE
The theory in the main text Fig. 4 is derived assuming a non-interacting Bose gas at thermal equilibrium in a symmetric, two-dimensional, harmonic trapping potential in the grand canonical ensemble [S8, S9]. We have also made use of an extension of the theory from spatial to temporal correlations [S10] which assumes that dissipative processes play no role. The theory is based on a series expansion in the fugacity. Fugacity is defined as ζ = exp (µ/k B T ), where µ is the chemical potential. The k th term in the expansion corresponds to occupancy of up to k particles in any given mode. Far below threshold, very few terms are needed for the series to converge. Arbitrarily large phase-space density can be obtained for negative chemical potential, i.e. ζ < 1, so the infinite series will always converge. The finite series also converges for large numbers of terms, albeit with a small positive chemical potential (ζ − 1 1) for large particle numbers. First-order correlations, normalised or otherwise, are calculated using equations (20)-(23) in Ref. [S10].
Coherence length is defined here as the size of a Gaussian function which fits the visibility as a function of shift V (r 0 , r 0 +xδ x , τ = 0), and not a fit to |g (1) (r 0 , r 0 +xδ x , τ = 0)| (the normalised first-order correlation function), which makes comparison to experimental data more robust against noise, especially in low-intensity regions of the images. Intensity scale is a Gaussian fit to the number density of photons as a S4 function of position. Number density can be extracted easily from the theory since number density at position r is equal to G (1) (r, r, τ = 0) (the un-normalised correlation function). Coherence time is the full-width at halfmaximum of g (1) (r 0 , r 0 , τ ). The correlation function in space matches well to a Gaussian, but the density does not. The temporal correlation function does not match any simple function (Gaussian, Lorentzian or symmetric double exponential decay), not least because, without damping, there are revivals of correlations functions at half-periods of the oscillations in the harmonic trapping potential. Real temporal coherence data is fitted with Gaussians.  S3. Demonstration that the series expansion used for the theory in the main text converges even above threshold. The numbers in the legend refer to the number of terms used in the series. In the main text, Fig. 4, 999 terms are used.
In Fig. S3 we show how well the finite series expansion converges, for coherence length, intensity length scale and coherence time. Far below threshold, the series converges as expected. Above threshold, the spatial scales converge when the series has a number of terms of the same order as the condensate population. The results qualitatively agree with exact calculations [S9].
The temporal correlation functions converge only below threshold. There is qualitative agreement with the dissipative model of Kirton and Keeling [S2, S11], in that coherence time does increase slightly with increasing number (comparable to the pumping rate) as threshold is approached. Since the model used takes no account of dissipation, it cannot predict anything but that the coherence time of a pure condensate ought to be infinite. The maximum coherence time inferred here is > 2.5 ps, limited by the range over which g (1) is evaluated, avoiding the non-observed revivials.
While there are no adjustable parameters in the theory, the theoretical photon number does not directly correspond to the experimental pump power. Below threshold, photon number and pump power are experimentally seen and theoretically expected [S11] to be linearly proportional. Likewise, far above threshold, but not so just above threshold. We can therefore trust our calculations only below threshold, and above threshold our calculations (assuming that the photon number remains proportional to the pump power) are plotted as dashed lines and are only a qualitative guide to what is expected. The calculations presented in the main text, Fig. 4, use 999 terms of the expansion.

A. Finite spatial resolution
The effect of finite imaging resolution and numerical aperture on measured interference patterns can be taken into account, starting from the known electric field at a point, E(r). Finite resolution is imposed by convolving with a point-spread function F (R): The overline indicates that finite resolution has been applied. The effect of finite numerical aperture is equivalent to applying the Fourier transform, applying a cutoff (multiplying by a top-hat function) and then inverse transforming, i.e. convolution with a cardinal sine, sin(x)/x. This function can then simply be absorbed in the definition of the point-spread function, F . The light at one output port of the interferometer is E P (r, r , τ ) = 1 √ 2 E(r, t) + E(r , t ) , where as usual τ = t − t . The effects of finite resolution are applied before the interference. The other output, Q, takes a minus instead of a plus. Then the intensity is: The first two terms are the intensities as seen with finite resolution. The last term written more explicitly is: It is possible to calculate the equilibrium first-order correlation function G (1) (r, r , τ ) for a non-condensed Bose gas in a harmonic trap if we make the strong approximations that there are no dissipative processes and that thermal equilibrium is respected [S10]. Using the symbols defined in Ref. [S10]: s is a co-ordinate x or y, K (k) s (s, s , t, t ) is the propagator, ζ the fugacity and β = 1/k B T the inverse temperature we obtain: where we have also assumed that point-spread function is separable: F (R) = F x (X)F y (Y ). This expression can be evaluated numerically, either by direct integration or via Fourier transforms. The finite-resolution correlation function is normalised: Calculated effect of finite imaging resolution on g (1) for the non-condensed photons, equivalent to P 5 × 10 −3 P th , using a 3 µm imaging resolution. The series expansion used 5 terms, which was sufficient for convergence so far below threshold.
The results are shown in Fig. S4 for a thermal cloud with pump power far below threshold (P 5×10 −3 P th ). The point correlation functions show shorter range coherence than those integrated over a finite resolution (a rotationally symmetric, 3 µm Gaussian point-spread function), and are consistent with the results seen at low pump powers in the main text, Figs. 2, 4 and 5.

S4. NON-EQUILIBRIUM THEORY OF PHOTON THERMALISATION AND CONDENSATION
Since the thermal-equilibrium theory clearly breaks down in the multimode regime, we have implemented the non-equilibrium model of Kirton and Keeling [S2, S11, S12]. To explain the multimode behaviour we are obliged to treat the spatial dependence of pumping. In Ref. [S12] it is suggested that the multimode behaviour might occur with inhomogeneous pumping, because of imperfect clamping of the excited-state fraction of dye molecules. To study the threshold for various modes, we evaluate Ref. [S12] equations (8) and (9). A crucial quantity in those equations is the absolute value of lightmatter coupling, denoted by Γ.
The light-molecule coupling strength can be approximately derived from measured quantities, noting that the typical timescale for population variations is: where N is the number of molecules, Γ dD and ρ dD are the light-matter coupling and molecular density in d dimensions and λ the wavelength of light in vacuum.
In three dimensions, the mean time between scattering events for photons moving at speed c * from molecules at volume density ρ 3D is τ typ = 1/ ρ 3D c * σ(λ). Here σ(λ) is the scattering crosssection at wavelength λ. Equating the timescales we find Γ max 3D = σ(534 nm)c * = 5.1 × 10 −12 m 3 /s. The variation of light-matter coupling with wavelength is known by interpolating experimental measurements enforcing a Kennard-Stepanov relations [S12]. In lower dimensions, the density is scaled typically by the cavity physical length and/or the harmonic oscillator length.

A. Multimode behaviour
The principal results are shown in the main text, Fig. 5, for the following parameters (symbols as used in Ref. [S12]) in one dimension: which is slightly smaller than the experiment (approx 10 l HO ), which compensates for the fact that we cannot efficiently perform the computations in two dimensions.
• Number of modes: 15, sufficient to make the main results converge.
• Cavity decay rate: κ = 10 9 s −1 • Molecular density: ρ 1D = 5.1 × 10 12 molecules/l HO which is equivalent to 1.7 mM solution concentration (similar to the experiments) if the extra dimension for conversion from 2D to 1D is taken to be the harmonic oscillator length, 6 µm. The length for conversion from 3D to 2D is the physical space between the mirrors, (q − q 0 )λ * /2 with q = 8 being the longitudinal mode number and q 0 4 expresses how far the electric field penetrates the surface of the dielectric mirrors.

B. Coherence length and time
The linewidth theory of Ref. [S11] takes into account only the single mode into which condensation occurs. Such a trunctation is not appropriate unless mode filtering optics are used. Instead, we make the approximation that dissipative processes are negligible, which means that we can write the classical electric field of the light in the cavity as a coherent sum over electric fields from the many modes: where the population of mode m is n m , its associated eigenfunction and angular frequency are ψ m (r) and ω m . The interferometer detectors can be used to measure the first order correlation function, much as in Eqn. (S6), but remembering that only stationary processes are observed. The result is: Note that since the electric field we wrote was classical, there is no quantum averaging here. Simple formulae that follow are: Results from these formulae are plotted in Fig. S5, summarised as the experimental data is, using Gaussian fits to extract coherence lengths and times, as well as realspace extents, G (1) (r, r, 0). For numerical reasons we cannot evaluate with enough modes to compare quantiatively to the experimental data, but the qualitative agreement is very good. For increasing pump powers, the spatial extent falls until threshold is reached, then rises in the multimode regime. Likewise the spatial coherence length rises then falls. The temporal coherence grows but above single-mode condensation threshold our approximation ignoring dissipative processes is not valid, until the multimode regime is reached. Agreement with the thermal-equilibrium model below threshold is good, and there is qualitative agreement with experimental results in all regimes (below, near and far-above threshold).