Topological order and equilibrium in a condensate of exciton-polaritons

We report the observation of the Berezinskii-Kosterlitz-Thouless transition for a 2D gas of exciton-polaritons, and through the joint measurement of the first-order coherence both in space and time we bring compelling evidence of a thermodynamic equilibrium phase transition in an otherwise open driven/dissipative system. This is made possible thanks to long polariton lifetimes in high-quality samples with small disorder and in a reservoir-free region far away from the excitation spot, that allow topological ordering to prevail. The observed quasi-ordered phase, characteristic for an equilibrium 2D bosonic gas, with a decay of coherence in both spatial and temporal domains with the same algebraic exponent, is reproduced with numerical solutions of stochastic dynamics, proving that the mechanism of pairing of the topological defects (vortices) is responsible for the transition to the algebraic order. Finally, measurements in the weak-coupling regime confirm that polariton condensates are fundamentally different from photon lasers and constitute genuine quantum degenerate macroscopic states.

tions on the basis of few and general parameters, the most important being dimensionality and symmetry [3][4][5] . The spontaneous symmetry breaking of Bose-Einstein condensates (BEC) below a critical temperature T > 0 is a remarkable example of such a transition, with the emergence of an extended coherence giving rise to a long range order (LRO) [6][7][8] . Notably, in infinite systems with dimensionality d ≤ 2, true LRO cannot be established at any finite temperature 9 . This is fundamentally due to the presence of low-energy long-wavelength thermal fluctuations (i.e. Goldstone modes) that prevail in d ≤ 2 geometries.
However, if we accept a lower degree of order, characterised by an algebraic decay of coherence, it is still possible to make a clear distinction between such a quasi-long-range-ordered (QLRO) and a disordered phase in which the coherence is lost in a much faster, exponential way. Such transition, in two dimensions (2D) and at a critical temperature T > 0, is explained in the Berezinskii-Kosterlitz-Thouless theory (BKT) by the proliferation of vortices-the fundamental topological defects-of opposite signs 10 . This theory is well established for 2D ensembles of cold atoms in thermodynamic equilibrium, where the transition is linked to the appearance of a linear relationship between the energy and the wavevector of the excitations in the quasi-ordered state 11 . The joint observation of spatial and temporal decay of coherence has never been observed in atomic systems, mainly because of technical difficulties in measuring long-time correlations. These are important observables to bring together because an algebraic decay, with the same exponent α, for both the temporal and spatial correlations of the condensed state, implies a linear dispersion for the elementary excitations [12][13][14] .
On the other hand, semiconductor systems such as microcavity polaritons (dressed photons with sizeable interactions mediated by the excitonic component) since the report of their condensation 15 appear to be ideal platforms to extend the investigation of many-body physics to the more general scenario of phase transitions in driven-dissipative systems 16 . However, establishing if the transition can actually be governed by the same BKT process as for equilibrium system has proven to be challenging from both the theoretical [17][18][19] and experimental perspective. Indeed, the dynamics of phase fluctuations is strongly modified by pumping and dissipation, and the direct measurement of their dispersion by photoluminescence and four-wave-mixing experiments is limited by the short polariton lifetime, by the pumping-induced noise and by the low resolution close to condensate energy [20][21][22] . Moreover, the algebraic decay of coherence has been experimentally demonstrated only in spatial correlations, while only exponential or Gaussian decays of temporal coherence, which are not compatible with a BKT transition, have been reported until now [23][24][25][26] . In particular, the lack of a power law decay of temporal correlations is a strong argument against a true BKT transition, as  Figure 1. Sketch of the pumping mechanism and interferogram. a, Spatial cross section (x,y) at different energies, where the vertical axes represents increasing energies. The carriers, injected by the pumping laser, relax quickly into excitonic states (yellow area) spatially confined within the pumping spot region.
Efficient scattering from the exciton reservoir into polariton states results in a region of high polariton density (red area) which expands radially outwards and relax into lower energy states. Above a threshold power, an extended 2D polariton condensate (blue area) is formed outside of the pumped region. b, Interference pattern of the extended 2D polariton condensate, obtained rotating one arm of the interferometer as the symbols on the right. The circles in yellow, dashed line indicate the blue-shifted region (emission is filtered around k=0) corresponding to the position of the laser spot.
will be demonstrated later on in a straightforward counter-example of a strongly out-of-equilibrium polariton system. For this reason, it has been a matter of constant debate in the field what is the nature of the various polariton phases, what are the observables that allow to determine a QLRO, if any, and how they compare with equilibrium 2D condensates and with lasers [27][28][29][30][31][32] . Recently, thanks to a new generation of samples with record polariton lifetimes, a thermalization across the condensation threshold has been reported via constrained fitting to Bose-Einstein distribution, suggesting a weaker effect of dissipation in these systems 33 . However, to unravel the mechanisms that drive the transition, and characterize its departure from the equlibrium condition, it is crucial to measure the correlations between distant points in space and time as we move from the disordered to the quasi-ordered regime 13,14,34,35 . So far, all attempts in this direction have been thwarted, not only by the polariton lifetime being much shorter than the thermalization time, but also by sample inhomogeneities 36,37 and/or small extents of the condensate. The latter, in particular, limits the power-law decay of the spatial coherence to the small spatial extension of the exciton reservoir set by the excitation spot 23,38,39 , which could result in an effective trapping mechanism 40 and finite size effects 28 .
In this work, using a high quality sample (in terms of long lifetimes and spatial homogeneity) to form and control a reservoir-free condensate of polaritons over a largely extending spatial region, we make the first observation (in any system) of the transition to a QLRO phase both in spatial and in temporal domains. Remarkably, the convergence of spatial and temporal decay of coherence allows us to identify the connection with the classic equilibrium BKT scenario, in which for systems with linear spectrum the exponents take exactly the same value α ≤ 1/4. 14 Stochastic simulations tuned to the experimental conditions, that reproduce the experimental observations in both space and time, further allow us to track vortices in each realisation of the condensate, confirming the topological origin of the transition. All these results settle the BKT nature of the 2D phasetransition for driven/dissipative polaritons in high quality samples, providing the equilibrium limit that allows to develop the paradigm of non-equilibrium regimes. Finally, this proves that a quasiordered polariton state is fundamentally different from a laser-excitons in a cavity under weak coupling conditions-for which a power-law decay of the first order coherence is observed only in space but not in time correlations, also putting to rest the long-time debate whether polaritons above threshold are distinct from a laser, and can indeed be called a condensate 41,42 . (blue data), stretched exponential (green data), power law (red data) spatial decay of g (1) (x, −x) and corresponding fitting residuals. d, e, f, g (1) (t, t + ∆t) for ∆t > 0 of the Gaussian (yellow data), stretched exponential (green data), power law (red data) temporal decay and corresponding fitting residuals. g, Blue line: β exponent from the stretched exponential fitting of g (1) (x, −x) as a function of polariton densities.
Red line: α exponent from the power law fitting of g (1) (x, −x) in the above threshold density range. h, Blue line: β exponent from the stretched exponential fitting of g (1) (t, t + ∆t) as a function of polariton densities.
Red line: α exponent from the power law fitting of g (1) (t, t + ∆t) in the above threshold density range. Error bars from fitting parameters errors.
The sample is excited non-resonantly (Methods), leading to the formation of the exciton reservoir (yellow region in Fig. 1a) which is localised within the pumping spot area due to the low exciton mobility. In turn, the repulsive interactions between excitons induce the energy blueshift of the polariton resonance at the center of the pumping spot (the dashed-white line in Fig. 1a is the energy blueshift in the x direction, corresponding to the Gaussian intensity profile of the exciting beam). Polaritons, which are formed beneath the exciton reservoir through energy relaxation, are much lighter particles than excitons and are accelerated outwards from the center of the spot. As g (1) g (1) 1.2 (green). Exponents α of the power law fitting for spatial (red) and temporal (orange). Note, we normalise the density d to the density d BKT at which the algebraic decay clearly fits the data better then the stretched exponential (see Fig. 6 in SI).
the 2D emission map with its centro-symmetric image. In Fig. 1b, the interference fringes are visible around the center of the image, which is the autocorrelation point r = r 0 . The phase coherence in space is therefore measured by the first order correlation function g (1) (r − r 0 ) (where omitted, r 0 = 0 and t = t is assumed). In Fig. 2, g (1) (r) is shown in colour-scale for different values of the polariton density. Increasing the pumping power, a higher level of coherence is sustained over larger distances ( Fig. 2(a-e)), up to a maximum in Fig. 2e. Note that the high quality of the sample allows a uniform coherence to extend over a wide spatial region of about 80 µm×60 µm. Increasing further the excitation power results in the shrinking of the spatial extension of the coherence, as shown in Fig. 2f.
In Fig. 3, the horizontal line profile |g (1) is studied at increasing pumping powers ( Fig. 3(a-c)). To allow a uniform description across the transition, both a power law and a stretched exponential functions are used in the fitting procedure: with B a scale parameter for the x-axis. At low pumping power (Fig. 3a), the decay is exponential and it is well fitted by eq. (2) with β ≈ 1. Increasing the pumping power, the exponential decay first changes into a short range envelope β > 1, and then, at the threshold density d th for the nonlinear increase of the ground-state population 43 , coherence starts to build up at larger distances. This corresponds to a change in the slope of the decay, which is still best fitted by eq. (2), but now with β < 1 (Fig. 3b). Finally, at d ≈ 2.7 d th , a power law decay is evident (Fig. 3c), with a high degree of spatial coherence (>50 %) extending over distances of ≈ 50 µm. Remarkably, this slow decay is characterized by the exponent α = 0.22. In Fig. 3g, the α and β exponents are shown at different densities: while below threshold only eq. (2) is able to fit the experimental data, above threshold the best fit is given by eq. (1) (see Supplementary Information). However, as will be shown in the following, it is essential to verify that a similar behavior is also observed for the temporal correlations.
The decay of temporal coherence is extracted by using the long delay line to change the relative temporal delay of the two interferometer arms up to 200 ps (see Methods and Supplementary Information). In Fig. 3(d-f), the temporal coherence at the autocorrelation point is shown for three different polariton densities. In Fig. 3h, the α and β exponents of equations (1) and (2) that best fit the experimental data are shown across the transition. Below threshold, coherence decays quickly and follows a Gaussian slope (β ≈ 2). At d = 1.3 d th , the temporal coherence can be best fitted by (2) with an exponent β ≈ 0.8 (or, with a slightly worst fit, with a power-law of exponent α ≈ 0.57), while at d ≈ 2.7 d th , the long time behaviour clearly follows a power law with α = 0.2. The residuals analysis proves the agreement between the experimental data and the fitting model (see Supplementary Information). Crucially, also for time correlations, α < 0.25, which coincides, within the experimental accuracy, with the one obtained from the spatial coherence at the corresponding density.
We performed complementary theoretical analysis, based on the exact solution of the stochastic equations of motions 19 (see Methods), with the same microscopic parameters as the ones of the experiment. We observed the same crossover from an exponential via stretched exponential to an algebraic decay of coherence in space and time (Fig. 4). In particular, for our long polariton lifetime, we see the spatial and temporal α being the same and always smaller then 1/4 above the BKT threshold (Fig. 4g), showing that the drive and dissipation do not prevail in this good quality sample in contrast to the earlier studied non-equilibrium cases 19,23 . Additionally, while the vortex-antivortex binding cannot be directly observed in the experiments, which average over many realizations, the numerical analysis is able to track the presence of free vortices in each single realisation, confirming the topological origin of the transition. Indeed, we see clearly that, in the algebraically ordered state, free vortices do not survive and the pairing is complete (Fig. 5 right). In contrast, the exponential and stretched-exponential regimes both show the presence of free vortices ( Fig. 5 left), the number of which decreases as we move across the transition. Since the stretched exponential phase is always associated with some presence of free vortices, this supports that we are observing a BKT crossover rather than a Kardar-Parisi-Zhang phase 17 . Finally, to highlight the importance of studying correlation not only in the space but also in the temporal domain, we used a microcavity in the weak-coupling regime, showing that a power law decay of coherence in space is associated in this case to a Gaussian decay in time. Using a sample with a lower quality factor and less quantum wells, we induce, under high non-resonant pumping, the photon-laser regime (as in a VCSEL) 44,45 . Despite the fact that this system is clearly out-of-equilibrium, it shows a power law decay of spatial coherence with α = 0.25 (Fig. 6a) within the pumping spot region (with a radius of about 10 µm). Remarkably, the behavior of spatial correlations is very similar to what obtained in Ref [23], but the temporal coherence, shown in

Conclusions
In conclusion, we observed the formation of an ordered phase in two-dimensional driven/dissipative ensemble of bosonic quasiparticles, exciton-polaritons in semiconductor microcavities, exploring both spatial and temporal correlations across the transition. This system lies at the interface between equilibrium and out-of-equilibrium phase transitions, and it has been often compared both to Bose-Einstein condensates and to photon lasers. We show that the measurement of spatial correlations g (1) (r) alone is not sufficient to establish whether an open/dissipative system is at equilibrium. Instead, two distinct measurements, one in time and one in space domain, are needed.
Satisfying this requirement, we report a power-law decay of coherence with the onset of the algebraic order at the same relative density and comparable exponents for both space and time correlations.
We should stress that the exceptionally long polariton lifetime in the present sample allows us to reach the BKT phase transition at lower densities, and in the region without the excitonic reservoir, resulting in a lower level of dephasing. In our experiments, the absence of any trapping mechanism, be it from the exciton reservoir or potential minima, allows us to avoid the influence of finite-size effects in the temporal dynamics of the autocorrelation 14 . Simulations with stochastic equations match perfectly the experimental results and demonstrate that the underlying mechanism of the transition is of the BKT type i.e. a topological ordering of free vortices into bound pairs, resulting in the coherence built up both in space and time. All these observations validate that polaritons undergo phase-transitions following the standard BKT picture, and fulfill the expected conditions of thermal equilibrium despite their driven/dissipative nature. Now that the equilibrium character of polaritons becomes a tunable parameter, the study of driven/dissipative phase transitions and of the universal scaling laws is within reach in this solid state device.

Methods
In the present experiment, polaritons were created in a planar semiconductor microcavity. The sample is described in detail in Ref. [43]  Gaussian spot with a FWHM ≈ 20 µm. The Gaussian laser profile in space is the best one to trigger the expulsion-relaxation mechanism, allowing the creation of a sufficiently uniform and extended polariton condensate. On the detection line, the emission from the sample is collected in epilayer configuration and it is sent to a Michelson interferometer with two delay lines as showed in Fig. S1 of SI. On the first interferometer arm, a short (fraction of the wavelength) piezo stage motorized is used to move a mirror. This allows the reconstruction of the spatial coherence g (1) (r − r 0 ) using the sinusoidal envelope of the intensity (see Supplementary Fig. S2). On the second arm, a long (about 200 ps) motorized screw actuator is used to change the position of the back-reflector to characterize the behavior of the coherence in time. By using the back-reflector, one of the two arms is centrosymmetrically rotated respect to the other and these two beams interfere on the focal plane of a CCD (Charge-Coupled Device) camera. To obtain energy resolved spatial and momentum emission images, a spectrometer is placed before the CCD detector. The coherence properties are investigated using a Michelson interferometer with two different delay lines, as showed in Fig. S1.

First-order Correlation Function Decay
The first order correlation function is defined as: with ψ * i and ψ i the creation and annihilation operators for the space-time point (r i , t i ), with i = 1, 2. By changing the length of one path through the use of a motorized piezo stage, the relative phase is changed, and the intensity of each point of the image shows a sinusoidal envelope as shown in Fig. S2. The path length is scanned by sub-wavelength steps in order to extract the first order correlation function as: where I ideal = (I 1 + I 2 )(2 I 1 I 2 ) −1 takes into account small asymmetries between the two interferometer arms, with I 1 and I 2 the intensities in the two interferometer arms and V is the visibility of the interference fringes obtained by fitting the data with: where the visibility V = A I 0 is shown in Fig. S2 and φ 0 is the initial phase.

Coherence Length and Time
In order to evaluate the coherence length and time as a function of density, both spatial g (1) (x, −x, t = 0) and temporal g (1) (t, t + ∆t) are fitted with stretched exponentials as reported in Eq. (2) of the main text. It is therefore possible to define both in space and time the relaxation length (time) as: with l e the renormalization factor scale of the x-axis and Γ the gamma function, with x-axis both the spatial and temporal one. The resulting coherence length and time are reported in Fig. S3:

Fitting Model and Residuals Analysis
The analysis of the fitting results allows to assess the applicability of the model used in the fitting procedure. In Fig. S4 we show the residuals of fitting the temporal first order correlation function g (1) (t, t + ∆t) in the quasi-ordered regime (Red square in Fig. 3h of the main paper). We report different attempts to fit the data using an exponential (blue, row (a)), a Gaussian (yellow, row (b)), a stretched exponential (green, row (c)) and a power law (red, row (d)). The P-P (Probability Probability) plot is a standard tool to investigate the deviance of a data set from the normal distribution.
Indeed, fitting residuals with a normal distribution around the zero value represents a strong indi- cation that the power-law model used to fit the data is the best one. The Kolmogorov-Smirnov test permits to quantitatively check the deviation from the normal distribution of a dataset. In order to assure that we have the minimum residual spreading, we used a normalised Gaussian distribution with σ the minimum from the residuals distributions for each model fitting the experimental data.
The value of this test, maximum for the power law model, combined with the fact (Fig. S4f) that the sum of the absolute values of the residuals is minimum for the power law model, confirms that the power law is the best fitting for the temporal decay of the correlation above BKT transition.
In Fig. S5 we show the same analysis for the same regime (red square in Fig. 3g of the main paper but for the spatial decay of correlations. In this case the Gaussian fitting is not shown because the values of the residuals is too large.

Theoretical Description
In the classic equilibrium BKT scenario, for a system with linear dispersion in the ordered phase, we expect a slow algebraic decay (as eq. (1) of the main text) of the first order coherence with exactly the same power-law exponents in both space (α s ) and time (α t ) S2,S3 , given by α s,t = k B T /n s , where T indicates temperature and n s the superfluid density. α s,t have an upper bound of 1/4 S4 at which density/temperature vortices proliferate, causing an exponential decay of coherence, characteristic for the disordered phase. At the same time non-equilibrium dissipative driven systems, with diffusive spectrum in the ordered phase S2,S5 , have been shown to still exhibit an algebraic decay of coherence but with temporal correlations decaying two times slower then the spatial ones α t = 1/2α s S2,S3 .
Moreover, values of α s as large as four times the equilibrium upper bound, when approaching the BKT transition, were reported both experimentally S6 and from theoretical analysis S19 , using beyond-mean-field truncated Wigner methods able to account for vortices, suggesting an "overshaken" superfluid state S7 . Finally, it has recently been suggested that the dissipation might in fact have an even more profound effect on the system with collective phase fluctuations destroying the algebraic order at long distances, leading to a stretched exponential decay of first order coherence characteristic of Kardar-Parisi-Zhang phase (KPZ) S8 . In that scenario the parameter β of the stretched exponential (see Eqn. (2) of the main paper) are also different for space (β ≈ 0.78) and time (β ≈ 0.48). Even if later estimates of the KPZ length-scales appeared to be unrealistic, and the presence of free topological defects strongly hampers the possibility of the KPZ phase S9 , the true nature of the 2D exciton-polariton phase transition and the resulting order is still at the center of an intense debate. Additionally, the type of the Renormalisation-Group (RG) analysis, which led to those conclusions S8,S9 , rely on the expansion in one over the mean-field density, and so are inadequate in describing the crossover region close to the phase transition. Thus, here, to specifically address the region close to the phase transition we restore to exact numerical solutions of the stochastic equations of motions, described in the next section.

Stochastic simulations
Our system consists of an ensemble of bosonic particles (the lower-polaritons) of mass m and lifetime 2κ, interacting via contact interactions g, and driven incoherently with pump of strength γ. The pumping saturation due to other processes is Γ. Using Keldysh field theory one can show that by including the classical fluctuations to all orders, but quantum fluctuations only to the second order (correct in the long-wavelength limit), and employing the MSR formalism, one can arrive at the stochastic equation for the field ψ(r, t) (for review see S10 ). Alternative derivation, using the Focker-Plank equation for the Wigner function, aimed at numerical implementation on a finite spatial grid dV has been reviewed in S11 . The finite grid version reads where dW is the Wiener noise with correlations where by |ψ(r, t)| − we abbreviated the following expression for the density |ψ(r, t)| − = |ψ(r, t)|− 1 dV , which comes from the Wigner function relating to the time symmetric and not time ordered operators. In Ref. S8 the Eqn. (S5) is solved approximately analytically using RG analysis by eliminating the amplitude fluctuations, and focusing solely on the non-topological phase-fluctuations. The procedure, however, relies on an expansion in one over the mean-field density, which together with approximating the amplitude fluctuations and discarding vortices, cannot capture the region close to the phase transition. Instead, here, we solve this equation exactly numerically on a finite grid.
In addition, we have checked our results with a more general model of a saturable drive, applicable There is no appreciable difference between the results from the two models for our relatively low densities.
We evolve the dynamics of the stochastic equations (S5) with the XMDS2 software S12 using a fixed-step (to ensure stochastic noise consistency) 4th order Runge-Kutta (RK) algorithm, which we have tested against fixed-step 9th order RK, and a semi-implicit fixed-step algorithm with 3 and 5 iterations. We choose the system parameters the same as in the experimental part: the mass of the microcavity lower polaritons is taken to be m C = 3.8 × 10 −5 m e , where m e is the electron mass, the polariton lifetime 2κ = 101.3ps, and the polariton-polariton interaction strength g = 0.004 meVµm 2 . The only parameter not possible to extract from the experiment is the saturation rate of the driving process (or in other words the three-body type losses in the system). We perform analysis for a range of Γ (or n s ) values, and since the other parameters are fixed, we choose the value of Γ (n s )to reproduce the overall length scale of g (1) (x).
In our method, the stochastic averages over the configurations of different realisations of the fields provide the expectation value of the corresponding symmetrically ordered operators, and it is important to get the results to converge in the number of realisations. The first order spatial correlation function g (1) (x) is evaluated according to Eq. (S1), by averaging over 100 independent stochastic paths, and additionally over auxiliary position in space r 0 , since in simulations our system is uniform. The temporal correlation function is evaluated from a single spatial point (to avoid picking up any spatial correlations) after the steady-state is reached, and averaged over 10000 stochastic paths. Since the Wigner average provides the expectation value of the corresponding symmetrically (and not time)-ordered operators, we need to subtract the expectation value of the commutator. For single time correlation functions, such as g (1) (x), the commutator is simply 1 2dV . For two-time correlation functions, such as g (1) (t), strictly speaking the expectation value of the commutator is unknown. It is, however, changing from 1 2dV at t = 0 to 0 at t → ∞. Using the two limiting values allows us to estimate the error, which we expect to be small given the densities considered. Indeed, for the densities used here the difference is practically indistinguishable. In 0.020 Residuals in , Figure S6. Residuals of the stretched exponential fit, β, for spatial (blue) and temporal (green) g (1) , and the power-law fit, α, for spatial (red) and temporal (orange) g (1) as a function of the polariton density normalised to the BKT threshold.
order to test the robustness of our conclusions to the choice of the numerical parameters, we perform simulations on different spatial grids varying from 2 to 5 µm in the grid spacing, and different system sizes from 256 to 1024 µm. Note, that to satisfy the condition necessary to derive the discrete version of equation (S5), g/[(κ + γ)dV ] 1, whilst maintaining a sufficient spatial resolution and, at the same time, a large enough momentum range, the window of available momentum grids is quite narrow. We find the main conclusions of our work (such as the crossover in g (1) (x) and g (1) (t) from exponential to power-law, spatial and temporal α being the same in the algebraic phase and smaller then 1/4, as well as the behaviour of vortices) independent on the choice of those parameters. Here, we present a case with 2 µm grid spacing, and system size of 256 µm, which is sufficiently larger from experimental to avoid any boundary effects influencing the relevant region. Finally, in order to assess which functional form fits the numerical data best we perform residuals analysis similar to those applied to experimental data. The residuals for the data and fittings presented in the main text is shown in Fig S6.