Non-equilibrium Berezinskii-Kosterlitz-Thouless Transition in a Driven Open Quantum System

The Berezinskii-Kosterlitz-Thouless mechanism, in which a phase transition is mediated by the proliferation of topological defects, governs the critical behaviour of a wide range of equilibrium two-dimensional systems with a continuous symmetry, ranging from superconducting thin films to two-dimensional Bose fluids, such as liquid helium and ultracold atoms. We show here that this phenomenon is not restricted to thermal equilibrium, rather it survives more generally in a dissipative highly non-equilibrium system driven into a steady-state. By considering a light-matter superfluid of polaritons, in the so-called optical parametric oscillator regime, we demonstrate that it indeed undergoes a vortex binding-unbinding phase transition. Yet, the exponent of the power-law decay of the first order correlation function in the (algebraically) ordered phase can exceed the equilibrium upper limit -- a surprising occurrence, which has also been observed in a recent experiment. Thus we demonstrate that the ordered phase is somehow more robust against the quantum fluctuations of driven systems than thermal ones in equilibrium.

The Hohenberg-Mermin-Wagner theorem prohibits spontaneous symmetry breaking of continuous symmetries and associated off-diagonal long-range order for systems with shortrange interactions at thermal equilibrium in two (or fewer) dimensions [2]. This is because long-range fluctuations due to the soft Goldstone mode are so strong as to be able to "shake apart" any possible long-ranged order. The Berezinskii-Kosterlitz-Thouless (BKT) mechanism (for an overview see Refs. [3,4]) provides a loophole to the Hohenberg-Mermin-Wagner theorem: Two-dimensional (2D) systems can still exhibit a phase transition between a quasi-long-range ordered phase below a critical temperature, where correlations decay algebraically and topological defects are bound together, and a disordered phase above such a temperature, where defects unbind and proliferate, causing exponential decay of correlations. Further, it can be shown [5] that the algebraic decay exponent in the ordered phase cannot exceed the upper bound value of 1/4.
The BKT transition is relevant for a wide class of systems, perhaps the most celebrated examples are those in the context of 2D superfluids, as in 4 He and ultracold atoms: Here, despite the absence of true long-range order, as well as a true condensate fraction, clear evidence of superfluid behaviour has been observed in the ordered phase [6]. Particularly interesting, and far from obvious, is the case of harmonically trapped ultracold atomic gases [7]. While, in an ideal gas, trapping modifies the density of states to allow Bose-Einstein condensation and a true condensate [8,9], weak interactions change the phase transition from normal-to-BEC to normalto-superfluid and recover the BKT physics despite the system's harmonic confinement.
These considerations are applicable to equilibrium, where the BKT transition can be understood in terms of the free energy being minimised in either the phase with free vortices or in the one with bound vortex-anti-vortex pairs. However, in recent years a new class of 2D quantum systems has emerged: strongly driven and highly dissipative interacting many-body light-matter systems such as, for example, polaritons in semiconductor microcavities [10], cold atoms in optical cavities [11] or cavity arrays [12,13]. Due to the dissipative nature of the photonic part, a strong drive is necessary to sustain a nonequilibrium steady-state. In spite of this, a transition from a normal to a superfluid phase in driven microcavity polaritons has been observed [14,15], and the superfluid properties of the ordered phase have started being explored [16][17][18][19][20]. Being strongly driven, the system does not obey the principle of free energy minimisation, and so it is not obvious whether the transition between the normal and superfluid phases, as the density of particles is increased, is of the BKT type i.e. due to vortex-antivortex pairs unbinding. Current experiments are not yet able to resolve single shot measurements, and so are not sufficiently sensitive to detect randomly moving vortices. Algebraic decay of correlations was reported from averaged data [1,21], however, the powerlaw decay displays a larger exponent than is possible in equilibrium, which posed questions as to the actual mechanism of the transition. On the theory side, by mapping the complex Ginzburg-Landau equation describing long-wavelength condensate dynamics onto the anisotropic Kardar-Parisi-Zhang (KPZ) equation, Altman et al. [22] concluded that although no algebraic order is possible in a truly infinite system, the KPZ length scale is certainly much larger than any reasonable system size in the case of microcavity polaritons.
In this work, we consider the case of microcavity polaritons coherently driven into the optical parametric oscillator regime [23,24] as the archetype of a 2D driven-dissipative system. Another popular pumping scheme is incoherent injec-arXiv:1412.7361v1 [cond-mat.quant-gas] 23 Dec 2014 tion of hot carriers, which relax down to the polariton ground state by exciton formation and interactions with the lattice phonons [14,15]. However, the incoherently pumped polariton system is challenging to model due to the complicated and not yet fully understood processes of pumping and relaxation. As a result, one is typically forced to use phenomenological models [25], which often suffer from spurious divergences. From this point of view, the parametric pumping scheme is particularly appealing, as an ab initio theoretical description can be developed in terms of a system Hamiltonian [26], and its predictions can be directly compared to experiment.
Analysing the non-equilibrium steady-state, we show that despite the presence of a strong drive and dissipation the transition from the normal to the superfluid phase in this lightmatter interacting system is of the BKT type i.e governed by binding and dissociation of vortex-antivortex pairs as a function of particle density, and bares a lot of similarities to the equilibrium counterpart. However, as recent experiments suggested [1], we find that larger exponents of the power-law decay are possible before vortices unbind and destroy the quasilong-range order leading to exponential decay. This suggests that the external drive, decay and associated noise favours excitations of collective excitations, the Goldstone phase modes, which lead to faster spatial decay, over unpaired vortices which would destroy the quasi-order all together. This externally over-shaken but not stirred quantum fluid constitutes an interesting new laboratory to explore non-equilibrium phases of matter.

Simulating driven-dissipative open systems
We describe the dynamics of polaritons in the OPO regime, including the effects of fluctuations, by starting from the system Hamiltonian for the coupled exciton and cavity photon field operatorsψ X,C (r, t), depending on time t and 2D spatial coordinates r = (x, y) ( = 1): Here, m X,C are the exciton and photon masses, g X the exciton-exciton interaction strength, and Ω R the Rabi splitting [10]. In order to introduce the effects of both an external drive (pump) and the incoherent decay, we add toĤ S a system-bath HamiltonianĤ SB whereψ l,k (t) are obtained Fourier transforming to momentum space k = (k x , k y ) the corresponding field operators in real spaceψ l (r, t).B l,k andB † l,k are the bath's bosonic annihilation and creation operators with momentum k and FIG. 1: Polariton system in the OPO regime. Upper panel: 2D map of the photonic OPO spectrum |ψ C,kx,ky =0 (ω)| 2 (logarithmic scale) of energy ω versus the kx momentum component (cut at ky = 0) for a single noise realisation and at a pump power fp = 1.02436f th p , where f th p is the mean-field OPO threshold. The arrows show schematically the parametric process scattering polaritons from the pump state into the signal and idler modes. Dashed (green) lines show the bare upper (UP) and lower polariton (LP) dispersions, while dotted (black) lines are the cavity photon (C) and exciton (X) dispersions. The solid (black) line underneath the spectrum is the ky = 0 cut of the single-shot time-averaged in the steadystate momentum distribution dt|ψ C,kx,ky =0 (t)| 2 , clearly showing the macroscopic occupation of the three OPO pump, signal, and idler states. Lower panels: 2D maps of the filtered space profiles |ψs,p,i(r, t)| 2 at a fixed time t for which a steady-state regime is reached -the pump emission intensity is rescaled to 1. Blue (red) dots indicate the vortex (antivortex) core positions.
energy ω l,k , describing the decay processes for both excitons and cavity photons. To compensate the decay, the system is driven by an external homogeneous coherent pump F (r, t) = f p e i(kp·r−ωpt) , which continuously injects polaritons into a pump state, with momentum k p and energy ω p .
Within the Markovian bath regime, standard quantum optical methods [27,28] can be used to eliminate the environment and obtain a description of the system dynamics in terms of a master equation. As the full quantum problem is, in practice, intractable, a simple yet useful description of the parametric oscillation process properties is provided by the mean-field approximation, where the quantum fieldsψ l (r, t) are replaced by the classical fields ψ l (r, t), whose dynamics is governed by the following generalised Gross-Pitaevskii equation [10] i∂ where κ X,C are the exciton and photon decay rates. By solving (1) both analytically and numerically, much work has been carried out on the mean-field dynamics for polaritons in the OPO regime and its properties analysed in detail [29][30][31]. Here, polaritons resonantly injected into the pump state, with momentum k p and energy ω p , undergo parametric scattering into the signal (k s , ω s ) and idler (k i , ω i ) statessee Fig. 1. As explained in detail in Ref. [32], as well as in other works [31], the full steady-state OPO photon emission ψ C (r, t) is filtered in momentum around the values of the signal, pump and idler momenta k s,p,i in order to get their corresponding steady-state profiles, i.e., ψ s,p,i (r, t) = |k−ks,p,i|<ks,p,i ψ C,k (t)e ik·r . The choice of each state filtering radius,k s,p,i , is such that the filtered profiles ψ s,p,i (r, t) are not affected by them -for details, see [32]. The meanfield onset of OPO is shown in the inset of Fig. 3, where the mean-field densities of both pump and signal are plotted as a function of the increasing pump power f p . At mean-field level, parametric processes lock the sum of the phases of signal and idler fields ψ s,i to that of the external pump, while allowing a global U (1) gauge symmetry for their phase difference to be spontaneously broken into the OPO phase -a feature which implies the appearance of a Goldstone mode [33]. As shown below, fluctuations above mean-field can lift, close to the OPO threshold, this perfect phase locking.
Fluctuations beyond the Gross-Pitaevskii mean-field description (1) can be included by making use of phase-space techniques -for a general introduction, see Ref. [34] , while for recent developments in quantum fluids of atoms and photons, see, e.g., Refs. [35][36][37]. Here, the quantum fieldsψ l (r, t) are represented as quasiprobability distribution functions in the functional space of C-number fields ψ l (r, t). Under suitable conditions, the Fokker-Planck partial differential equation, which governs the time evolution of the quasiprobability distribution, can be mapped on a stochastic partial differential equation, which in turn can be numerically simulated on a finite N × N grid with spacing a (along both x and y directions) and a total size L x,y = N a comparable to the polariton pump spot size in state-of-the-art experiments. For the system under consideration here, the Wigner representationone of the many possible quasiprobability distributions -is the most suitable to numerical implementation: in the limit g X /(κ X,C dV ) 1, where dV = a 2 is the cell area, it appears in fact legitimate [26,38] to truncate the Fokker-Planck equation, retaining the second-order derivative term only, thus obtaining the following stochastic differential equation: (2) Here, dW l=X,C are complex valued, zero-mean, independent Weiner noise terms with dW * l (r, t)dW m (r , t) = δ r,r δ l,m dt dV , and the operator H M F coincides with H M F in Eq. (1) with the replacement |ψ X | 2 → |ψ X | 2 − 1 dV . The same stochastic equation (2) can be alternatively derived applying a Keldysh path integral formalism to the Hamiltonian H S +Ĥ SB , integrating out the bath fields, and keeping only the renormalisation group relevant terms [39]. Note that, remarkably, some of the difficulties of the truncated Wigner method met in the context of equilibrium systems, such as for cold atoms, are suppressed here by the presence of loss and pump terms, i.e., the existence of a small parameter g X /(κ l dV ) which controls the truncation [10]. Note, however, that the bound on this truncation parameter involves the cell area dV of the numerical grid, that is the UV cut-off of the stochastic truncated Wigner equation. For typical OPO parameters, it is possible to choose dV small enough to capture the physics, but at the same time large enough to keep the UV issues under control.
We reconstruct the steady-state Wigner distributions ψ l (r, t) by considering a monochromatic homogeneous continuous-wave pump F (r, t) = f p e i(kp·r−ωpt) as before and letting the system evolve to its steady-state. In order to rule out any dependence on the chosen initial conditions, we have considered four extremely different cases: empty cavity with random noise initial conditions and adiabatic increase of the external pump power strength; mean-field condensate initial conditions; either random or mean-field initial conditions in the presence of an unpumped region at the edges of the numerical box, so to model a sort of vortex-antivortex reservoir. The different initial stage dynamics, and their physical interpretation for each of these four different initial conditions, are carefully described in [32]; in all four cases we always reach the very same steady-state regime, i.e., all noise averaged observable quantities discussed in the following lead to the same result -this could not be a priori assumed for a non-linear system.
Below we first analyse results from single noise realisations (concretely, here, for the case of mean-field initial conditions and no "V-AV reservoir" present), by filtering the photon emission at the signal, pump and idler momenta as also previously done at mean-field level. The filtered profiles are again indicated as ψ s,p,i (r, t) -for details on filtering see [32]. Second, we consider a large number of independent noise realisations and perform stochastic averages of appropriate field functions in order to determine the expectation values of the corresponding symmetrically ordered quantum operators. In particular, we evaluate the signal first-order correlation func-tion as (3) where the averaging . . . is taken over both noise realisations as well as the auxiliary position R, and where t is either a fixed time after a steady-state is reached, or we take additional time average in the steady-state [32].

Vortices and densities across the transition
It is particularly revealing to explore the steady-state profiles, i.e. ψ s,p,i (r, t), of the signal, pump and idler states. Fig. 1 shows a cut at k y = 0 of the OPO spectrum, |ψ C,kx,ky=0 (ω)| 2 , determined by solving Eq. (2) for ψ X,C (r, t) to a steady-state and evaluating the Fourier transforms in both space and time. Note, that the logarithmic scale of this 2D map plot (which we employ to clearly characterise all three OPO states) makes the emission artificially broad in energy, while in reality this is sharp (as required by a steadystate regime), as well as it is very narrow in momentum. The filtered space profiles ψ s,p,i (r, t) shown in the bottom panels of Fig. 1 reveal that while the pump state is homogeneous and free from defects, vortex-antivortex (V-AV) pairs are present for both signal and idler states. Note, that while at the meanfield level the sum of the signal and idler phases is locked to the one of the pump (and thus a V in the signal implies the presence of an AV at the same position in the idler), the large fluctuations occurring in the vicinity of the OPO threshold make this coherent phase-locking mechanism only weakly enforced, resulting in a different number (and different core locations) of V-AV pairs in the signal and idler states. Because the density of photons in the idler state is much lower than the one at the signal (see, e.g., the photonic momentum distribution plotted as a solid black line inside the upper panel of Fig. 1), while both states experience the same noise strength, the number of V-AV pairs in the filtered photonic signal profile is much lower than the number of pairs in the filtered photonic idler profile. Phase locking between signal and idler is recovered instead for pump powers well above the OPO threshold, where long-range coherence over the entire pumping region is re-established.
The proliferation of vortices below the OPO transition, followed by a sharp decrease in their density and their binding into close vortex-antivortex pairs is illustrated in Fig. 2. Here, we plot the 2D maps of the phase for the single noise realisation of the filtered OPO signal ψ s (r, t) (photonic component) for increasing values of the pump power f p in a narrow region close to the mean-field OPO threshold f th p ; the position of the generated vortices (antivortices) are marked with blue (black) dots. While at lower pump powers there is a dense "plasma" of Vs and AVs, the number of V-AV pairs decrease with increasing pump powers till eventually disappearing altogether (not shown). We do also record a net decrease in the distance between nearest neighbouring vortices with opposite winding number with respect to that between vortices with the same winding number. In order to quantify the vortex binding across the OPO transition, we measure, for each detected vortex, the distance to its nearest vortex, r V-V and to its nearest antivortex r V-AV ; and similarly, for each detected antivortex, we measure r AV-AV and r AV-V . We then consider the symmetrised ratio b = rV-V+rAV-AV rV-AV+rAV-V . In order to extract a noise realisation independent quantity, an average over many different realisations, as well as over individual vortex positions, is performed to obtain b ; this quantity b → 1 for an unbound vortex plasma, while b → 0 when vortices form tightly bound pairs. We observe a dramatic drop in b (green squares in Fig. 3) when increasing the pump power across the OPO threshold, indicating that vortices and antivortices are indeed binding, as it is expected for a BKT transition.
By evaluating other relevant noise averaged observable quantities, we are able to construct a phase diagram for the OPO transition in Fig. 3 and link it with the properties of the BKT transition. We evaluate the averaged signal photonic density at some time t in the steady-state, n s = dr |ψ s (r, t)| 2 /V , where V = (N a) 2 is the system area and . . . indicates the noise average for the stochastic dynamics (blue dots). We also show the steady state signal density in the mean field (orange squares). The corresponding mean-field densities for both signal (orange line) and pump (black line) are presented for comparison in the inset of Fig. 3. At mean-field level, both signal and idler (not shown) suddenly switch on at the OPO threshold pump power, f p = f th p and both states are macroscopically occupied above threshold. The effect of fluctuations is to smoothen the sharp meanfield transition, as clearly shown by the (blue) dots in the main panel of Fig. 3, where we plot the noise average signal density n s . This is because, even below the mean-field threshold, incoherent fluctuations weakly populate the signal. Note also that, even though somewhat smoothened, we can still appreciate a kink in the n s density, but at higher values of the pump power compared to the mean-field threshold f th p . We identify this as the novel BKT transition for our out-of-equilibrium system, as discussed more in detail below. This is further confirmed by a sudden decrease of the averaged number of vortices in the signal (red diamonds), and of the averaged distance between nearest neighbouring vortices of opposite winding number, b (green squares), as a function of the pump power concomitant with the observed kink for n s . These results suggest that the system undergoes an OPO transition which, by including fluctuations above mean-field, is indeed analogous to the equilibrium BKT transition. Both vortices and antivortices proliferate below some threshold and, above, they bind to eventually disappear altogether. As indicated by the black square in the inset of Fig. 3, the region for such a crossover is indeed narrow in the pump strength.

First-order spatial correlations
For systems in thermal equilibrium, the BKT transition is associated with the onset of quasi-off-diagonal long-range order, i.e., with the algebraic decay of the first-order correlation function in the ordered phase, where vortices are bound, and exponential decay in the disordered phase, where free vortices do proliferate. In order to investigate whether the same physics applies to our out-of-equilibrium open-dissipative polariton system, we evaluate the signal first-order correlation function g (1) (r) according to the prescription of Eq. (3) and characterise its long-range behaviour in Fig. 4. We observe the ordering transition as a crossover in the long-distance behaviour between an exponential decay in the disordered phase, g (1) (r) ∼ e −r/ξ , and an algebraic decay in the quasi-ordered phase, g (1) (r) ∼ (r/r 0 ) −α . We therefore fit the tail of the calculated correlation function to both of these functional forms and observe that, at the onset of vortex binding-unbinding and proliferation, the signal's spatial correlation function changes its long-range nature, from exponential at lower pump powers to algebraic at higher (see Fig. 4).
However, in contrast with the thermal equilibrium case, we do observe that the exponent α of the power-law decay (inset of Fig. 4) can exceed the equilibrium upper bound of 1/4 [5], and can reach values as high as α 1.2 for f p /f th p = 1.0136, just within the ordered phase. Further, as thoroughly discussed in Ref. [32], it is interesting to note that, close to the transition, we do observe a critical slowing down of the dynamics: Here, the convergence to a steady-state is dramatically slowed down compared to cases above or below the OPO transition, a common feature of other phase transitions. At the same time, close to threshold, the convergence of noise averaged number of vortices is much faster than the convergence of the power-law exponent α. This indicates that the fluctuations induced by both the external drive and the decay preferentially excite collective excitations, rather than topological excitations, resulting in vortices being much more dynamically stable. Finally, note that for sufficiently strong pump powers, the power-law exponent becomes extremely small and thus quasi-long-range order is difficult to distinguish from the true long-range order over the entire system size.
Our findings explain why recent experimental studies, both in the OPO regime [40], as well as for non-resonant pumping [1], experienced noticeable difficulties in investigating the power-law decay of the first order correlation function across the transition. We do indeed find that the pump strength interval over which power-law decay can be clearly observed is extremely small, and the system quickly enters a regime where coherence extends over the entire system size, as measured in [40]. This was also observed in non-resonantly pumped experiments when using a single-mode laser [14]. However, by intentionally adding extra fluctuations with a multimode laser pump, as in [1], power-law decay was finally observed in the correlated regime, with an exponent in the range α 0.9 − 1.2, in agreement with our results.

Discussion
Using microcavity polaritons in the optical parametric oscillator regime as the prototype of a driven-dissipative system, we have numerically shown that a mechanism analogous to the BKT transition, which governs the equilibrium continuous-symmetry-breaking phase transitions in two dimensions, occurs out of equilibrium for a driven-dissipative system of experimentally realistic size. Notwithstanding the novelty and significance of this result, there are a number of novel features which warrant discussion as they are peculiar to non-equilibrium phase transitions. We have observed that the exponent of algebraic decay in the quasi-long-range ordered phase exceeds what would be attainable in equilibrium. This recovers a recent observation [1], and strongly suggests that indeed a non-equilibrium BKT may have been seen there. Moreover, our findings imply that the ordered phase is more robust to fluctuations induced by the external drive and decay than an analogous equilibrium ordered phase would be to thermal fluctuations. Although for realistic experimental conditions, we have found that the region for BKT physics, before the pump power is strong enough to induce perfect spatial coherence over the entire system size, is indeed narrow, we believe our work will encourage further experimental investigations in the direction of studying the non-equilibrium BKT phenomena. Even though the small size of the critical region has so far hindered its direct experimental study, our calculations indicate that the macroscopic coherence observed in past polariton experiments [10,14,15,23,24,40,41] results from a non-equilibrium phase transition of the BKT rather than the BEC kind.

Methods
We simulate the dynamics of the stochastic equations (2) with the XMDS2 software framework [42] using a fixed-step (where the fixed step-size ensures 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 to be close to current experiments [18]: The Rabi frequency is chosen as Ω R = 4.4 meV, the mass of the microcavity photons is taken to be m C = 2.3 × 10 −5 m e , where m e is the electron mass, the mass of the excitons is much greater than this so we may take m −1 X → 0, the exciton and photon decay rates as κ X = κ C = 0.1 meV, and the exciton-exciton interaction strength g X = 0.002 meVµm 2 [43]. The pump momentum k p = (k p , 0), with k p = 1.6 µm −1 , is fixed just above the inflection point of the LP dispersion, and its frequency, ω p − ω X (0) = 1.0 meV, just below the bare LP dispersion. In order to satisfy the condition necessary to derive the truncated Wigner equation (2), g X /(κ X,C dV ) 1, whilst maintaining a sufficient spatial resolution and, at the same time, a large enough momentum range so that to resolve the idler state, simulations are performed on a 2D finite grid of N × N = 280 × 280 points and lattice spacing a = 0.866 µm. Thus, the only system parameter left free to be varied is the pump strength f p : We first solve the mean-field dynamics (1) in order to determine the pump threshold f th p for the onset of OPO. We then vary the value of f p around f th p in presence of the noise in order to investigate the nature of the OPO transition. We analyse the results from single noise realisations by filtering the full photonic emission for signal pump and idler, as described in the main text, as well as in [32]. Further, we average all of our results over many independent realisations, which are either taken from 96 independent stochastic paths or from multiple independent snapshots in time after the steady-state is reached: As thoroughly discussed in [32], 96 stochastic paths is shown to be sufficient to ensure the convergence of noise averaged observable quantities. For each noise realisation, vortices are counted by summing the phase difference (modulo Supplemental Material for "Non-equilibrium Berezinskii-Kosterlitz-Thouless Transition in a Driven Open Quantum System" In this supplementary information we provide a more detailed account of some of the technical issues arising from our numerical methods, which are not of sufficiently general interest to warrant a discussion in the main text, but which are germane to the validity of our conclusions. In particular, we discuss the filtering we perform in order to isolate the signal state, and the loss in resolution this incurs, and we describe various tests we have undertaken in order to ensure that the results we present are properly converged.

I. MOMENTUM FILTERING
In order to study the condensate at the signal (or idler), we have to filter the full emission ψ C (r, t) in such a way as to omit contributions with momentum outside a set radius about the signal (idler) states in k−space. This applies both to the theory and the experimental data, and this procedure could equivalently be performed in either momentum or energy space. Here we filter in momentum space, and define: where k s,p,i is the momentum of the signal, pump and idler states respectively andk s,p,i is the filtering radius. We fix the filtering radius for the signal and idler to bek s,i = 1 2 |k s,i −k p | and choose to operate in the frame of reference co-moving with the signal (idler) so that the (physically irrelevant) background current does not show in the data. The filtering reveals the phase-freedom of the signal (idler) state at the expense of spatial resolution. Note, that this limits our ability to distinguish vortex-antivortex (V-AV) pairs with a separation less than a distance π/k l , which for our chosen parameters is ≈ 2.708µm. However, such extremely close vortex-antivortex pairs do not affect the spatial correlations of the field at large distances.

A. Number of stochastic realisations
The stochastic averages over the configurations of different realisations of the fields provide the expectation value of the corresponding symmetrically ordered operators. The realisation-averaged results presented in the main paper have been averaged over 96 realisations taken from independent stochastic paths as we assessed this to be sufficiently many to give a reliably converged result. In Fig. 1 we show the first order correlation function g (1) (x) for f p = 1.017f th p , averaged over different numbers of realisations (96,192,288,384 and 480). We conclude that 96 is sufficient to determine the nature of the correlations.
The average number of vortices also does not change significantly (no more than ±0.5) beyond 96 realisations, and so we conclude that this number is sufficient to ensure consistent results in both the smooth and topological sectors of the model.

B. Convergence in time to a steady state
In order to assess when a time evolution has reached its steady state we consider the average (over realisations) of both the signal density, |ψ s | 2 , and the number of vortices. In Figs. 2 and 3 we show the evolution of these quantities toward a steady state, assuming the initial condition wherein the system is allowed to reach its mean field steady state (at time t = 0) and stochastic processes are adiabatically switched on. In practice this means that the noise terms are multiplied by a ramp function such that the noise is slowly increases from 0, achieving half its maximum at t 0 , at a rate determined by t r . For the empty cavity with noise initial condition, in which the pump is adiabatically increased, it is the pump term that is multiplied by this ramp function. We choose t 0 = 450ps and t r = 150ps as these prove to give a fairly rapid convergence to the steady state. It is, however, important to note that the steady state is unique and does not depend upon the specific values of t 0 and t r we choose nor on the initial conditions, which we start the dynamics from. In all cases the steady state is reached within around twenty nanoseconds. The slowest convergence occurs in the vicinity of the BKT transition. This critical slowing of the dynamics is to be anticipated given the divergence of the correlation length as the transition is approached. C. Realisations from independent stochastic paths vs. independent time snapshots in a single path In Fig. 4 we show the number of vortices for a broad range of pump powers across the transition averaged over realisations taken from independent stochastic paths (red dots) and from multiple snapshots over time (blue dots) once the steadystate is reached. Averaging over time within the steady state is a less numerically intensive approach, and shows an excellent agreement with the average over realisations even in the critical region. In practice, averaging over time is computationally efficient away from the critical region, where the steady state is reached quickly and memory is quickly lost during the time evolution. Around the critical region it takes a significant time to reach a steady state and then to decorrelate the snapshots, so averaging over stochastic paths becomes more computationally effective. The blue dots show the number of vortices averaged over the final 500 frames, running from 15ns to 30ns, of a single stochastic path, once the system has converged to a steady state. We observe an excellent agreement between these two approaches and so away from the critical region we can consider the (less numerically intensive) average over time.

D. Dependence on different initial conditions
We have investigated a number of physically diverse initial conditions in order to rule out dependence of the final steady state on the starting configuration. To initiate our simulations we either adiabatically increase the pump strength atop a white noise background or adiabatically increase the strength of the stochastic terms starting from the mean field steady state. We also wish to eliminate the possible effects of trapped vortices, and so in addition to simulations with a uniform pump region we have performed simulations, where a small strip around the numerical integration box is left unpumped so as to act as a source/sink for vortices. The four distinct initial conditions we consider are therefore those classified in the following In schemes A and D, there are initially very many vortices, which then proceed to annihilate with one another so that the overall vorticity tends to decrease toward the steady state. In schemes B and C there are no vortices at the outset but as stochastic processes shake the system the vorticity increases toward the steady state. Therefore even before the true steady state is reached, we can treat schemes A and D as upper bounds for the vorticity and schemes B and C as lower bounds. In schemes C and D we incorporate a reservoir (region with no drive and thus of very low density) of width 34.72nm along the two sides of the numerical integration box parallel to k p , which we then exclude from calculations of the vorticity, signal density, and correlations. The width of the reservoir does not change the result except in that it needs to be wide enough that the condensate has room to decay away to zero.
In Supplementary Movie 1 we show the evolution towards the steady state for all four of these schemes with a pump power f p = 1.017f th p . We plot the vortex density (top panel) as well as dynamics of vortices (bottom panel, where blue and red dots show positions of V and AV cores) for each scheme (A, B, C and D starting from the left). We see that every scheme converges fairly rapidly towards the same steady state.
From the animations it is evident that the vortices always appear in pairs but that these pairs are not always tightly bound. It was not a priori obvious that all four schemes should lead to the same steady state but we take this as evidence that the physical process leading to the steady state is universal for this system. We observe that the schemes incorporating an empty reservoir converge more quickly to the steady state than their counterparts without a reservoir. Nevertheless we present data for scheme B because the reservoirs have no counterpart in the traditional BKT transition to which we wish to compare our results.

III. FITTING TO THE SPATIAL CORRELATION FUNCTION
The behaviour of the first order correlation function g (1) (x), whether it decays as an exponential or as a power law, is determined by fitting the tail of the data to each functional form and taking the better fit of the two. The fitting is illustrated in Fig. 5. For pump powers close but above the critical region it is clear that the power-law fits the data better (which is also confirmed by a larger correlation coefficient), while below the critical region the exponential fit is closer to the data. For stronger pump powers, away from the transition, g (1) (x) is practically constant on these length-scales (as reported in experiments and discussed in the main text).