Dynamical critical exponents in driven-dissipative quantum systems

We study the phase-ordering of parametrically and incoherently driven microcavity polaritons after an infinitely rapid quench across the critical region. We confirm that the system, despite its driven-dissipative nature, fulfils dynamical scaling hypothesis for both driving schemes by exhibiting self-similar patterns for the two-point correlator at late times of the phase ordering. We show that polaritons are characterised by the dynamical critical exponent z ~ 2 with topological defects playing a fundamental role in the dynamics, giving logarithmic corrections both to the power-law decay of the number of vortices and to the associated growth of the characteristic length-scale.

Two-dimensional (2d) quantum driven-dissipative systems, such as exciton-polaritons in microcavities [1], display rich universal critical phenomena due to the interplay between drive, dissipation and potential anisotropy and their collective dynamics induced by the Bosedegeneracy [2].The picture is already quite complex even at the level of the steady-state, where the system has been predicted to adopt Kardar-Parisi-Zhang (KPZ) [3] or Berezinskii-Kosterlitz-Thouless (BKT) [4] type scaling, depending on the subtle interplay between drive and dissipation, and finite size and effective anisotropy [5,6].Moreover, the even more advanced question of how continuous drive and dissipation affect dynamical critical behaviour remains largely open.
In this vein, one of the most useful concepts is the scaling hypothesis.Complex systems at criticality display self-similar patterns, since the system is statistically equivalent after a global and arbitrary change of scale [7].The scaling phenomenon also reveals the existence of universal critical exponents, characterising macroscopic properties of the system at critical points.Such critical exponents for driven-dissipative microcavity polaritons have not been measured to date while theoretical predictions are subject to debate.Specifically, 2d drivendissipative systems described by KPZ-like phase dynamics are expected to show the dynamical critical exponent z ≈ 1.61 while in the equilibrium limit, where the KPZ nonlinearity ceases to be important, we expect z = 2.Moreover, it has been suggested that the dynamical critical exponent z for microcavity polaritons coupled to a reservoir takes values of either 1 or 2, depending on system parameters [8], putting in question the universality of the phase ordering in this system.
In this letter, we explore the critical properties of 2dmicrocavity polaritons, and in particular the character-istic length and density of topological defects (vortices), by studying the phase ordering (scaling) dynamics after an infinitely rapid quench from the disordered to deep in the (quasi)ordered [7,9] phase.In order to capture the universal properties of the phase ordering process we study the polariton system under three different pumping configurations, focusing on experimentally realistic conditions.Specifically we study the coherent pumping in the optical parametric oscillator (OPO) regime, and the incoherent pumping (IP) with and without a frequencyselective pumping mechanism.In particular, we report strong numerical evidence that the three different polariton systems all show the same critical behaviour, characterised by the dynamical critical exponent z ≈ 2 with logarithmic correction to the diffusive dynamics.Thus, we find that the universal properties of the system, also in the dynamical case, are dominated by BKT-type of physics in analogy to the static case considered recently theoretically [10] and experimentally [11].It seems that the KPZ nonlinear terms do not play a role during dynamical crossing of the critical point, at least for the realistic system sizes considered in this work.
System and method.
Our system consists of an ensemble of bosonic particles (excitons (X) and photons (C) for parametrically-pumped, and lower polaritons (LP) for incoherently-driven case) with finite lifetime interacting via contact interactions in two dimensions, and driven in two distinct ways: parametrically and incoherently.The dynamical equations can be derived using Keldysh field theory by including the classical fluctuations to all orders, but quantum fluctuations to the second order, which is appropriate in the long-wavelength limit, and employing the Martin-Siggia-Rose (MSR) formalism (for review see [2]).An alternative derivation can be performed using Fokker-Planck equations for the Wigner function trun-cated to the third-order [1,12].Both methods lead to the same stochastic equation for the field ψ(r, t) with the noise term accounting for quantum and thermal fluctuations.
For the parametrically driven case, the finite grid version with = 1 reads [10]: where ψ X,C = ψ X,C (r, t) (with r = (x, y)) are the exciton, and cavity-photon fields respectively, dW X and dW C are the complex-valued zero-mean white Wiener noise terms with dW * l (r, t)dW m (r , t) = δ r,r δ l,m dt.The external monochromatic coherent pump F p = f p e i(kp•r−ωpt) injects photons with momentum k p and frequency ω p while both fields decay with their corresponding rates κ X and κ C : here we use a shorthand notation , with m X and m C (= 2.3 10 −5 m e ) the exciton and cavityphoton masses respectively.Since m X m C , we consider the limit m X → ∞ and, consequently, the exciton field kinetic energy term disappears from Eq. ( 1).The exciton-photon Rabi-splitting is given by Ω R , the exciton-exciton interaction strength by g X and dV = a 2 is the numerical grid unit cell area (lattice spacing a).We solve Eq. ( 1) for parameters typical of current OPO experiments [1,10,13], namely: Ω R ≈ 4.4meV , g X ≈ 2 × 10 −3 meV µm 2 .We consider κ X = κ C , with κ C = 1/τ C and photon lifetime τ C = 6.58ps.We set k p = (1.6,0)µm −1 , ω p = ω LP (k p ) where ω LP is the lower polariton dispersion.Here, we present results for a grid of 512 2 lattice points with lattice spacing a = 0.87µm (total unit cell of L x = L y = 444.42µm).
Since we are interested in the universal properties of driven-dissipative systems, we also consider the alternative typical set-up of incoherent driving.Under the assumption that the high-energy reservoir follows adiabatically the condensate evolution and that the exciton and photon are locked into a single lower-polariton branch, the equation reads ( = 1) [14] with ψ LP = ψ LP (r, t) and the lower-polariton field density reads (after subtracting the Wigner com-  n s .We use typical experimental parameters [15]: γ LP = 1/τ LP with the polariton lifetime τ LP = 4.5ps, polariton mass m LP = 6.2 10 −5 m e , polariton-polariton interaction strength g LP = 6.82 10 −3 meV µm 2 .To improve the physical relevance of the model when approaching the critical region, following [14,16] we implement frequency-selective pumping mechanism so that relaxation to low-energy modes is favoured over energies higher than the cut-off frequency ω cut Ω/ 1 + |ψ LP | 2 − /n s which are now suppressed.We use here P/P th = 1.06,where P th = γ LP is the mean field critical pump, Ω = 50γ LP = 11.09ps−1 and n s = 500µm −2 (labelled as IP Ω=50 ), for the frequency-dependent pumping, and P/P th = 1.1, Ω = ∞ and n s = 1500µm −2 for the frequency-independent driving (labelled IP Ω=∞ ).We solve Eq. ( 2) in a square lattice of 301 2 points, lengths L x = L y = 295.11µmand periodic boundary condition and average over a sufficientlylarge number (400) of realisations in all schemes.The healing lengths ξ = 1/ 2m ξ fulfilled in all cases [17].In both models a non-equilibrium steady state with finite particle density |ψ X,C,LP | 2 is established once the pumping strength overcomes the cavity losses.Also by tuning the strength of the pump power both the OPO and the IP system (with/without frequency-dependent pumping mechanism) undergo a BKT-type of phase transition between a disordered phase with exponential and an ordered phase with power-law decay of spatial coherence, governed by the binding/unbinding of vortexantivortex pairs (see Fig. 1).For the parametric pumping, the BKT transition in the steady-state is analysed in detail in [10].A similar transition takes place for IP system (Fig. 1 bottom).It is worth noting that in this case, stronger fluctuations at higher modes (Ω → ∞) and smaller saturation density (n s → 0) lead to a larger shift of the pump power of the BKT transition with respect to the mean-field onset of macroscopic population growth.
Universal dynamical scaling.We now study an infinitely-rapid quench across a critical point.For the parametrically pumped case (see Fig. 1 top), we quench through the upper critical threshold f up p .The system is prepared in the steady-state of a deeply disordered phase at a given pump power f i p > f up p where bound and unbound vortices proliferate (see the inset), and is instantaneously quenched to a deep quasi-ordered regime by adjusting the external drive to f f p , with f f p < f up p .For the incoherently pumped case (see Fig. 1 bottom), we quench from the random initial configuration deep in the disordered phase (i.e.steady-state for P i = 0), to a quasi-ordered regime by setting the external drive to P f > P th at t = 0 and letting the system evolve.We explore the dynamical scaling properties of the system during the phase ordering by considering the first order two-point correlation function [10]: where . . .denotes averaging over both noise realisations and the auxiliary position u, and t is the time after the quench.For the parametrically pumped scheme, ψ in Eq. ( 3) corresponds to the signal, which is obtained by filtering the cavity-photon field ψ C from Eq. ( 1) around the signal momentum k s .For incoherently pumped system, ψ is the polariton field ψ LP .

SS
IP Ω = 50 0 50 100 150 /g (1) SS )(r) First-order correlation function g (1) (normalised by the corresponding steady state correlator g SS ) for the parametricallypumped (top) and incoherently-pumped polaritons with (bottom) frequency-dependent pumping as a function of the rescaled distance r/L(t) at different times during the phaseordering process.The apparent collapse of the curves confirms the dynamic scaling hypothesis.Insets show g (1) /g (1) SS at different times, from which the characteristic length-scale L(t) is obtained by considering (g (1) /g SS )(L(t), t) = 0.5 (white dots).
In the phase ordering kinetics of the planar XYmodel in 2-dimensions [18][19][20], the non-steady state twopoint correlator (3) fulfils the dynamical scaling form: SS (r, t) • F (r/L(t)), with the steady state correlation function decaying algebraically at long distances, as g SS (r) ∼ r −α .The scaling function F tends to 1 when r L(t), indicating that the critical correlations have been established at distances much smaller than L(t) at time t, which defines the characteristic lengthscale of the system L(t).Since our system is highly nonequilibrium it is far from obvious whether similar scaling behaviour holds here in the presence of strong drive and dissipation.Indeed, we obtain a perfect collapse when plotting the two point correlation function g (1) /g SS as a function of the rescaled length r/L(t) at different times of the late dynamics for both driving schemes (see Fig. 2).For consistency, we extract the length-scale L(t) when (g (1) /g SS )(L(t), t) = 0.5 (white dots in insets of Fig. 2), with the independence of our conclusions on the intersection value verified in the Supplementary Material [17].Since our polariton system fulfils the scaling hypothesis, as revealed by the collapse of the two-point correlation function, we can access the universal dynamical critical exponent z of our driven-dissipative system by analysing the growth of the characteristic length L(t) and the decay of the number of topological defects (vortices) at late times after a sudden quench.Note, that an equilibrium analogue of the phase degree of freedom for polariton system is the planar XY-model, where free vortices and bound vortex-antivortex pairs exist even in the steady-state below and above the BKT phase transition respectively.The existence of the steady-state vortices plays a fundamental role in the phase ordering process and introduces the characteristic logarithmic correction into the late-time dynamics both of L(t) and number of vortices n v following a sudden quench [18,21], such that L(t) ∼ ((t/t 0 )/ log(t/t 0 )) 1/z and n v (t) ∼ ((t/t 0 )/ log(t/t 0 )) −(2/z) where t 0 is a nonuniversal microscopic system timescale (taken here as t 0 = 1ps [17]).The two relations follow from the fact that n v (t) ∼ 1/L(t) 2 when there is a unique length-scale in the system, which is true at late times in the dynamics [17].
To demonstrate the emergence of this universal scaling from our numerical data, Fig. 3 plots the characteristic length L(t) and number of vortices n v (t) for the three different pumping schemes considered.OPO (top) and IP Ω=∞ (centre) cases are qualitatively identical in the sense that g (1) collapses in the same time window as that in which the respective vortex dynamics reaches the converged z ≈ 2 value.We are careful to fit L(t) and n v late enough in the quench so the dynamics becomes universal (indeed both reveal z ≈ 2) but before size effects, power-law correlations or very small vortex number affect our analysis.However, in the case of IP Ω=50 (bottom), strong damping of collective fluctuations, introduced by the explicitly frequency-dependent nature of the pump, is responsible for the collective modes to reach z ≈ 2 at a much earlier time than the topological modes.Nevertheless, both channels show the dynamical critical exponents to be z ≈ 2 in their appropriate late time windows.A more detailed discussion about the fitting criteria can be found in Supplementary Material [17].
We stress that the sufficiently late-time analysis is essential to allow all channels to equilibrate properly, and fulfil the scaling hypothesis, whereas fitting early in the phase-ordering process, and before the dynamics becomes universal (light green regions in Fig. 3), can lead to the incorrect conclusion of z ≈ 1.Note that for the experimentally realistic parameters considered in our simulations, the phase-ordering takes place on timescales similar to the polariton population growth [17], which is the case for which Ref. [8] predicted z ≈ 1.  Summary and outlook.We have extended the study of universal critical properties, specifically the dynamical critical exponent, of strongly driven-dissipative twodimensional quantum systems by considering an infinitely rapid quench from the disordered to a quasiordered phase.Our work reveals clearly that such universal properties of the system only emerge for appropri-ately late time dynamics of both the topological defects and smooth phase fluctuations, whereas the early dynamics are nonuniversal and can give misleading information about the value of the critical exponents.
Importantly, for all pumping configurations we find a value z ≈ 2 for the dynamical critical exponent, different from the one of conservative Bose systems [22][23][24][25].On the other hand, our results indicate that under realistic polariton pumping schemes and experimental parameters, the relevant scaling is analogous to the one predicted for the planar XY-model [18][19][20] and, therefore, the non-equilibrium physics brought by the KPZ nonlinearity does not show in the phase ordering of realistic size systems.
Supplementary Material for: Dynamical critical exponents in driven-dissipative quantum systems In this Supplementary Material we present technical details related to our numerical procedure and analysis, offering conclusive proof that z ≈ 2 with logarithmic corrections is the correct dynamical exponent for excitonpolariton systems across all experimentally-relevant regimes.After explicitly demonstrating numerical convergence, and the different physical regimes probed in terms of the interplay of density-growth and vortex-decay dynamics, we demonstrate the clear emergence of z ≈ 2 in the presence of logarithmic corrections at sufficiently late evolution times; we also show the direct relation between vortex number and obtained characteristic length scale, demonstrate the effective independence of our conclusions on the choice of intersection point for the correlation function collapse and show that careful consideration of the conventionally-ignored intrinsic nonuniversal time-scale t 0 does not affect our findings.

NUMERICAL METHODS AND CONVERGENCE
First order correlation function.For both parametric and incoherently pumped case, the first order correlation function g (1) (r) is calculated by means of Eq. ( 3) in the main text.
For the incoherently pumped case, we first compute the two-point correlation function for each vertical and horizontal array of a single realisation of the wave function ψ.The resulting N × N outcomes are then averaged and normalized before additionally averaging over N p realizations, and computing corresponding error bars.For the parametric pump case, we filter the full emission ψ C , in such a way as to omit contributions with momentum outside a set radius about the signal states, in momentum space with a sharp step-like rectangular filter; for more details see [S1].The noise-averaged autocorrelation function of the filtered field ψ s is computed efficiently in momentum space.
Vortex Number Calculations.The number of vortices is evaluated from the phase gradients around closed paths of each grid point.We extract the phase, φ, and gradient of the phase, ∇φ = v, using finite differences.A grid point is identified as having a vortex when the circulation Γ = C v dr 2π.
For the incoherently pumped system, a Gaussian (lowpass) filter is first applied to the wave function, in order to remove all the high frequency noise components.This removes all noise with wavelength (in pixels) smaller than, or of the order of, the standard deviation of the filter's Gaussian kernel.
Convergence in lattice size.Fig. S1 provides conclusive evidence of the independence of our results on lattice size for both OPO and IP regimes, based on which the analysis in our main paper is conducted for grid spacings a = 0.87µm (OPO), a = 0.98µm (IP) and grid sizes 444.42 µm (OPO) and l = 295.11µm(IP).Specifically, the top (middle) panel show a comparison of the timeevolution of the vortex density for the OPO (IP) case upon increasing the total grid size, while keeping dis-cretisation fixed.Fig. S1 (bottom) shows corresponding IP results (for Ω = 50γ LP ) for the correlation function, with the coloured bands indicating the error bars; this demonstrates that our numerics is sufficient to avoid the correlation function being affected by boundary effects even as the system approaches the critical region.
Convergence in stochastic numerical realisations.Having fixed the grid spacing and size, we next investigate the dependence of results on realisation over different numerical trajectories (done in addition to averaging over the entire grid in single realisations).Although Fig. S2 shows that on first inspection a number of ≈ 100 (OPO) or 50 (IP) realisations may be enough for sufficient convergence in the overall dynamics, we nonetheless stress that averaging over a high number of realisations is essential for an accurate determination of the long-time dynamical critical exponent for L and n v .Throughout our work we have therefore used N p = 400 realisations for both OPO and IP cases.
Convergence in computational method (Incoherently pumped system).To ensure independence of our numerical results on computational method, we have implemented two distinct schemes for the incoherentlypumped case, based on finite differences method for the Laplacian implemented in Fortran and based on spectral methods and the publicly-available XMDS2 code [S2].Fig. S3 (top) shows excellent agreement between the two cases even for only 50 numerical realisations.Throughout our work, we have chosen to report results based on the XMDS2 spectral method, as this has exponential convergence, i.e. the error scales as ∝ a N with increasing resolution (number of grid points N) -as opposed to the algebraic scaling ∝ a 2 of finite difference schemes.
Convergence in cut-off choice (Incoherently pumped system).Within the stochastic Gross-Pitaevskii equation, the grid spacing a sets a maximum momentum cutoff ∝ a −1 in the numerical representation of the quantum field, restricting the maximum number of evolving modes.It is therefore important to check that our results do not depend on the chosen cut-off value.Fig. S3 (bot- tom) demonstrates that for two different grid spacings a = 0.49 µm and a = 0.98µm (while keeping the total box size fixed to l = 295.11µmby varying the number of grid points) the evolution of n v is practically identical.

VORTEX NUMBER DYNAMICS
A previous study [S3] of incoherently-pumped excitonpolariton systems has observed the dynamical critical exponent to be dependent on the quality of the sample, through the polariton lifetime, thus arguing for different types of non-universal dynamics for the driven-dissipative exciton-polariton system.Contrary to such a statement, our work demonstrates unequivocally that the universality class of the phase transition in polariton systems of FIG.S3: Convergence in methods and cut-off choice for the incoherently pumped system.Top panel: evolution of number of vortices for two different computational schemes: spectral methods (red points) and finite differences methods (blue points).Bottom panel: time evolution of number of vortices for two different lattice spacings, i.e. cutoff choice, a = 0.49µm (red points) and a = 0.98µm (blue points).Parameters as in Fig. S1.
experimentally realistic size falls, as anticipated, within the 2D-XY-model universality class provided that exponents are extracted at the appropriate late times.
Fig. S4 displays time evolutions of the averaged particle density and average number of vortices for the OPO (top) and IP (middle/bottom) systems.During the growth of the density of the degenerate exciton-polariton system, the vortex pairs decay monotonically in time, eventually reaching the long-time limit where universal phase-ordering kinetics features should set in.
Ref. [S3] has effectively argued that the value of the dynamical critical exponent depends on the interplay between time-scales for density saturation and vortex pair annihilation, indicating that in cases where the two processes occur effectively in parallel one would expect the critical exponent z = 1, characteristic of conservative superfluids [S4, S5], whereas in the case where the majority of vortex pair annihilation occurs after the density has effectively saturated, one would expect z = 2, as in the 2D XY model [S6, S7].
Fig. S4 shows the three cases considered within our numerics, based on experimentally-relevant parameters [S8, S9].Specifically we consider cases in which the two IP Ω = ∞ FIG.S4: Evolution of averaged density and number of vortices.Density growth and pair annihilation for coherent (top), frequency-dependent (central) and frequencyindependent incoherent pumping systems (bottom).We report regions where vortex dynamics follows Eq. ( S3) with critical exponent z ≈ 1 (faint green region) and z ≈ 2 (grey region) as reported in Fig. 3 main paper.Parameters as in main paper.
processes occur in parallel (bottom), the density growth is initially much faster than vortex annihilation (middle), and the density grows after most vortex annihilation has taken place (top).As reported in the main text (Fig. 3), when fitting all such cases with the applicable scaling law with logarithmic corrections we find clearly that z ≈ 1 at early evolution times (highlighted in Fig. S4 by the faint green regions), and that the dynamical critical exponent always converges to z ≈ 2 at later times (grey regions in same plots), as convincingly shown in next section.
Here t 0 is a nonuniversal microscopic system time-scale and consideration of all relevant system time scales (see last section) demonstrates our results to be insensitive to its precise value.
For completeness, Fig. S5 depicts a characteristic ex-FIG.S5: Phase ordering dynamics of an IP system.Snapshots capturing the phase ordering process for condensate density in µm −2 (top) and phase (bottom) of an incoherently pumped system.The initially noisy configuration of the system (left) becomes gradually ordered through the annihilation of vortex pairs (centre and right) according to the scaling law of Eq. (S5).Parameters as in Fig. S1.
ample of density and phase profiles during the coarsening process.

LOGARITHMIC CORRECTIONS AND DYNAMICAL EXPONENT
We now discuss the fitting criteria we adopt in the latetime stage dynamics (dark grey regions in Fig. 3 in main paper) of the coarsening process.The time evolution of the average number of vortices during the annealing process is fitted with different curves; • A power law formula for a diffusive system • a function with logarithmic corrections to account for the dynamics of pairs of topological defects: where A 0 , z and z are free parameters of the fit, and z, z correspond to the extracted dynamical critical exponents in the presence, or absence, of logarithmic corrections.
Constraining the dynamical exponent to the anticipated z = 2 value [S7], we also show fits to defects dynamics with Comparison of different fitting curves are shown in Fig. S6 for OPO (top), IP Ω=∞ (centre) and IP Ω=50 (bottom) cases within our region of convergence (grey band), from which we can infer numerous conclusions.Firstly, while to lowest order all fits of Eqs. ( S2)-(S3) provide a reasonable description, careful consideration rules out the value z = 2 in the absence of logarithmic corrections (black dashed lines).Moreover, although fits of Eqs.(S2) and (S3) (solid green and yellow lines respectively) are practically indistinguishable (even in terms of their residuals!),we stress that the exponent z extracted by fitting Eq. (S3) [with logarithmic corrections] is always closer to the anticipated theoretical value for the dynamical exponent (z = 2) than the extracted exponent z obtained from the fit based on Eq. (S2).This suggests that, to the extent that the dynamical exponent should be consistent with a value of 2, logarithmic corrections have to be present in the system, thus confirming the theoretically-anticipated description [S6, S7].A de- IP Ω = ∞ , n v y ∼ (t/t 0 ) −2/z y ∼ ((t/t 0 )/log(t/t 0 )) −2/z FIG.S8: Same as Fig. S7 but for frequency-independent pumping scheme.L(t) from our simulations by intersecting the graph of g (1) /g (1) SS (L(t)) at a value of η = 0.5, ensuring we are insensitive to both short-range and finite-size effects.We then fit the extracted L(t) by the corresponding relations [S6, S7] From such plots, corresponding to different parameter space and pumping schemes, we can infer the following general conclusion about the evolution of the exponent z: • Early time dynamics with z ≈ 1 2: At early times, the system is rapidly approaching (and eventually crossing) the critical region.The dynamical exponents from both characteristic length L and vortex number n v agree, but scaling hypothesis does not hold in this regime (i.e. the correlation function do not collapse onto each other).Within this early-stage evolution, it is possible to identify a precise time-window where z ∼ 1 (shown by the vertical light green bands in Figs.S7-S9).
• Intermediate dynamics with z ≈ 2 but no discernible logarithmic corrections: During this stage, the system has crossed the transition but not relaxed enough for logarithmic corrections to become visible.Both functions are well-fitted by Eqs. ( S2) and (S7), yielding a dynamical exponent z ∼ 2 which do not include logarithmic correction.
• Late-time dynamics with z ≈ 2 with evident presence of logarithmic corrections: After sufficient evolution, the system enters the regime where the scaling hypothesis holds: in this regime the logarithmic corrections are measurable and the dynamical exponent z ∼ 2 is consistent with fitting curves with Eqs. ( S3) and (S8).At such times, the vortex decay has slowed down, with the system left with vortex-antivortex pairs only.
We stress that in all cases, the exponents extracted at late times by fittings curves with logarithmic corrections are much closer to the expected value of z = 2 than corresponding fits without logarithmic corrections.As a guide to the eye, Figs.S7-S9 highlight the interval 1.9 < z < 2.1 by a horizontal light-blue shaded region.
Summarising, we note that our detailed analysis demonstrates that late-time dynamics of parametrically and incoherently pumped exciton-polariton systems follow a unique scaling law, consistent with that theoretically predicted for two-dimensional geometries in the context of the 2D XY model [S6, S7].We also stress that achieving the correct dynamical exponent z = 2 requires very extensive numerical simulations featuring long-time evolution (such that the system fully enters the correlation-function collapse window where the scaling hypothesis holds), a very high temporal resolution during all dynamics and a large number of independent numerical realisations.

SELF-CONSISTENCY OF CHARACTERISTIC LENGTH-SCALE
The scaling hypothesis implies the existence of a unique length scale L(t).This length scale extracted from the collapsing correlation functions should agree with the corresponding length scale L v (t) obtained from the vortex density through L v (t) ∼ n v (t) −1/2 which coincides to the mean distance between vortices.Such a correspondence should happen at sufficiently late times in the system evolution.
Fig. S10 demonstrates clearly that, at appropriately late times the scaled quantity ñv L2 , where ñv = n v (t)/n v (t max ) and L = L(t)/L(t max ), converges to a constant value of one, confirming that the mean dis- Analysis on the dependence of zL on the intersection point η.We study the dependence of exponents z extracted from the characteristic length scale L(t) on intersection condition (g (1) /g SS )(L) = η within the collapse region for the three systems analysed: OPO (red), IPΩ=50 (blue) and IPΩ=∞ (green).Values of exponents zL (solid lines) are always closer to 2 (grey band) than exponents zL (dashed lines).
tance between vortices L v (t) is always proportional to the growing length scale L(t) during the dynamic scaling.

DEPENDENCE OF zL ON THE INTERSECTION POINT η
In order to obtain the characteristic length-scale from the correlation function essential for demonstrating dynamical scaling, we have used the condition (g (1) /g (1) SS )(L) = η, with the value of η chosen as 0.5 (OPO, IP Ω=50 and IP Ω=∞ ).
Here we confirm the validity of our results, by a detailed consideration of the behaviour of the exponents z L and zL extracted by fitting L(t) (with curves Eqs.(S7) and (S8) respectively) over the entire time window for which the scaling hypothesis holds (corresponding to the dark grey regions of Fig. 3 in main paper) as a function of the intersection point η.Such a dependence is shown in Fig. S11, which confirms our key point that the exponent z L obtained in the presence of logarithmic corrections is closer to the predicted value of 2, than the numerically-extracted exponent, zL , obtained in the absence of logarithmic corrections.
Moreover, we stress that the exponent z L lies within the range 1.9 < z L < 2.1 for a wide range of η values including values of η used to extract the growing length scale L(t) in previous works, i.e. η = 0.25 in Ref. [S3] and η = 0.3 in Ref. [S7].Deviations only occur at quite extreme values of η, where finite-size and short-range effects come into play.The introduction of the anticipated logarithmic corrections through the related formulas n v (t) = A 0 • [(t/t 0 )/log(t/t 0 )] −2/z L(t) = B 0 • [(t/t 0 )/log(t/t 0 )] 1/z (S9) introduces a nonuniversal microscopic system timescale, t 0 .The presence of t 0 is typically ignored, on the grounds that it should not affect universal features, but its detailed consideration is nonetheless relevant for the precise slope of the corresponding curves for n v and L. Following common practice [S18, S19], throughout this work we have effectively chosen t 0 such that it "drops out" of the corresponding equations.By considering the dependence of the universal dynamical exponent z on this nonuniversal parameter, we effectively shed more light onto why such an approach is acceptable.Since t 0 is a typical system timescale, it should be of broadly the same order of magnitude as all other relevant system timescales.These are listed, for both OPO and IP schemes in 10 −1 ps < τ system < 10 1 ps.This already suggests that the conventional choice of t 0 taking the value of 1 in the system units is acceptable, which for our current purposes would imply the already-chosen value of t 0 = 1ps.However, in order to be certain that our results do not critically depend on our choice of t 0 , Fig. S12 displays the dependence of the dynamical exponent z (in the presence of logarithmic corrections) on t 0 within the anticipated time window 10 −1 ps < t 0 < 10 1 ps.From this figure it is evident that, within this range, all our numericallyextracted values for z take a value 1.9 < z < 2.1 (depicted by the horizontal light grey area in Fig. S12).We thus conclude that our extensive numerical analysis confirms that all three physical cases considered (OPO, IP with and without frequency-dependent pumping) are consistent with the presence of logarithmic corrections and a dynamical critical exponent of z = 2, as expected for the system to be in the same universality class as the 2D XY model.

FIG. 1 :
FIG.1: BKT transition with parametric and incoherent pumping.Top panel: noise-averaged density from stochastic equations of the signal (blue), idler (red) and pump (green) as a function of pump power fp for parametrically driven polaritons across OPO.Bottom panel: meanfield (dotted lines) and noise-averaged (solid lines) densities for frequency-independent (green) and frequency-dependent (blue) incoherent pumping.In both panels arrows indicate the infinitely rapid quench protocol across the critical region and pump powers are scaled to their corresponding mean-field threshold values.The insets show typical snapshots of the real space phase profile for the initial and late-time states.

FIG. S1 :
FIG. S1: Convergence in lattice size.Top panel: convergence in time of vortex density for OPO system with fp/f up p = 0.97, for two different box of size l = 222.208µm and l = 444.42µm with respectively N = 256 and N = 512 number of points and grid spacing a = 0.87 µm.Central panel: convergence in time of vortex density for IP system with P/P th = 1.5, Ω = 11.09ps −1 and ns = 500 µm −2 , for two different box of size l = 295.11µm,l = 444.42µm with respectively N = 601 and N = 901 number of points.Bottom panel: spatial convergence of first order correlation function for IP system before it enters the power-law stage at t = 4.5 ns for the two different boxes.

FIG. S2 :
FIG.S2: Convergence in number of stochastic paths.Top panel: convergence in time of number of vortices for OPO system averaged over Np = 100, 200, 400 stochastic paths.Central panel: convergence in time of number of vortices for IP system averaged over Np = 50, 400 noise realizations.Bottom panel: spatial convergence of first order correlation function for IP system before it enters the power-law stage at t = 4.5 ns over Np = 100, 200, 400 stochastic paths.Parameters as in Fig.S1.
IP Ω = ∞ , L y ∼ (t/t 0 ) 1/z y ∼ ((t/t 0 )/log(t/t 0 )) FIG.S9: Same as Fig.S7but for frequency-dependent pumping scheme.Note the different time axis in this case associated with the disparity in time scales for collective and topological channel equilibration.See Fig.3(bottom) in main text.
) where B 0 , z and z are free parameters.The time evolution of the numerically-extracted dynamical exponents are shown in Figs.S7-S9 for both cases of vortex (top) and length scale (bottom) determination, with extracted average exponents indicated by circles; the widths and heights of the "bands" around such points correspond respectively to the temporal averaging intervals and the error bars from the numerical fits.The vertical green and grey bands indicate corresponding regions shown in Figs. 3 of the main paper.
FIG. S10:Comparison of the number of vortices with the characteristic length scale of the system.Scaled number of topological defects ñv = nv(t)/nv(tmax) times the scaled characteristic length L = L(t)/L(tmax) squared as a function of time for the parametrically (blue) and incoherently pumped (red) system.We show results from the start of the quench (t = 0) to the end of the collapse region (t = tmax).We observe that at late-time dynamics the system fulfils ñv L2 ∼ 1, indicating that during dynamical scaling the mean domain size Lv ∼ n −1/2 v is proportional to the length scale L(t).
FIG.S12: Dependence of z on the time-scale t0.Dependence of exponent z obtained by fitting both number of vortices and characteristic length within the collapse region for the three cases.