Kardar-Parisi-Zhang universality in discrete two-dimensional driven-dissipative exciton polariton condensates

The statistics of the fluctuations of quantum many-body systems are highly revealing of their nature. In driven-dissipative systems displaying macroscopic quantum coherence, as exciton polariton condensates under incoherent pumping, the phase dynamics can be mapped to the stochastic Kardar-Parisi-Zhang (KPZ) equation. However, in two dimensions (2D), it was theoretically argued that the KPZ regime may be hindered by the presence of vortices, and a non-equilibrium BKT behavior was reported close to condensation threshold. We demonstrate here that, when a discretized 2D polariton system is considered, universal KPZ properties can emerge. We support our analysis by extensive numerical simulations of the discrete stochastic generalized Gross-Pitaevskii equation. We show that the first-order correlation function of the condensate exhibits stretched exponential behaviors in space and time with critical exponents characteristic of the 2D KPZ universality class, and find that the related scaling function accurately matches the KPZ theoretical one, stemming from functional Renormalization Group. We also obtain the distribution of the phase fluctuations and find that it is non-Gaussian, as expected for a KPZ stochastic process.


I. INTRODUCTION
Scale invariance and universality are the most striking features unveiled by the study of equilibrium continuous phase transitions.They refer to the fact that in the vicinity of a critical point, systems of completely different nature can exhibit correlations in their order parameter fluctuations characterized by the same power-law, with the same universal critical exponents.
The emergence of universality is well understood in systems at thermal equilibrium, but far less so out of equilibrium.In such systems, criticality can, for instance, emerge spontaneously during the time evolution of the system, a phenomenon known as self-organized criticality [1].A paradigmatic example of this feature is the mechanism of kinetic roughening of a stochastically growing interface [2,3], for which the celebrated Kardar-Parisi-Zhang (KPZ) equation was originally derived [4].It was realized later on that a strikingly large variety of systems actually belong to the KPZ universality class [5], such as switching fronts propagating in liquid crystals [6], or the development of modern urban skylines under building constraints [7] to cite a few.
Recently, KPZ universality was found in a surprising context: in the phase dynamics of driven-dissipative exciton polariton (polariton) condensates.Polaritons are nonequilibrium quasi-particles that form within semiconductor microcavities operating in the strong coupling regime between bound electron-hole pairs (excitons) and cavity photons.Under non-resonant excitation, polaritons can spontaneously condense into a macroscopic coherent state [8][9][10].In this state, the microscopic model describing the dynamics of the condensate phase has been shown to map onto the KPZ equation [11][12][13].This has far reaching conceptual implications, as it shows that the nature of a polariton condensate is very different from its thermal equilibrium counterpart.KPZ universality was confirmed numerically in a one-dimensional (1D) polariton condensate, both in terms of scaling of spatiotemporal correlations [14], and in terms of the phase fluctuation statistics [15,16].In a polariton condensate, in contrast with more usual KPZ interfaces, space-time vortices can emerge stochastically due to the intrinsic 2π periodicity of the phase [17].However, KPZ dynamics can co-exist with -and remain robust to -such defects.Indeed, KPZ scaling in a 1D polariton condensate has been reported experimentally [18], thus promoting polaritons as a compelling experimental platform to probe KPZ universality.
In two dimensions (2D), far fewer results have been obtained for the KPZ equation: no exact result is available, and mainly functional Renormalization Group calculations [19][20][21][22] and numerical simulations [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39] have provided quantitative information.Experimentally, polaritons offer a unique opportunity to investigate KPZ physics, as they are well suited to form a two-dimensional condensate.The dynamics of its phase in 2D enables to study the 2D KPZ universality class releasing the stringent requirement of realizing a growing surface embedded in a three-dimensional physical space.In this work we address theoretically and numerically the emergence of KPZ physics in a non-resonantly driven discrete 2D polariton condensate.Our work is complementary to the studies on coherent excitations [40][41][42], the optical parametric oscillator [43], and the quadratic driving scheme [44], where the possibility of observing KPZ-like behavior was highlighted.All of the above correspond to different models with respect to the one considered here.
In 2D, equilibrium bosonic quantum fluids undergo the Berezinskii-Kosterlitz-Thouless (BKT) transition stemming from vortex-antivortex unbinding.Owing to the driven-dissipative nature of polaritons, the behavior of vortices is different: the attractive vortex-anti-vortex interaction is, for instance, dramatically damped at large distances, and can even change sign to become repulsive [45,46], hindering their annihilation.In some specific non-equilibrium conditions, the vortex motion can also be self-accelerated, leading to the creation of additional vortex pairs [45,46].The role of vortex unbinding in driven-dissipative condensates have also been studied using the duality with a theory for nonlinear electrodynamics coupled to charges [47,48].Thus, vortex-filled phases are also present in non-equilibrium systems and could be detrimental to the emergence of KPZ universality in 2D.This is supported by an argument based on perturbative Renormalization Group (RG), which suggests that the proliferation of vortices, which happens at spatial scale L v , spoils KPZ stretched exponential decay of coherence characterized by a much longer length scale [13].Moreover, several studies, based on both experiments and numerical simulations of 2D polariton condensates, have reported BKT-like types of behaviors, with specificities that are related to the non-equilibrium regime [49][50][51][52][53].
An important feature of this system is that the dimensionality of the parameter space is large, such that a multitude of realistic regimes dominated by KPZ could still exist.For example, the momentum-dependent loss rate of the polaritons plays an important role in the magnitude of the KPZ non-linearity λ, and hence in the possibility of reaching a KPZ regime.Indeed, KPZ scaling of spatiotemporal correlation functions were then found in some region of the parameter space, as well as vortexfree regimes [54].Yet, this was obtained in the low noise limit, while a realistic noise input is typically larger and fixed by the polariton pump and losses.
In this work, we demonstrate the existence of a KPZ regime for the phase dynamics of polariton condensates in 2D, using, for most parameters, values that can be realistically achieved in experiments.We show that in this regime, the proliferation of topological defects such as spatial or spatiotemporal vortices can be limited by acting on the pump strength and by introducing a cutoff length scale via discretization of space thanks to the introduction of a lattice potential.Polariton lattices are routinely created in semiconductor microcavities (see eg [55] for a review).We obtain the KPZ scaling in the spatiotemporal correlation function of the condensate wavefunction.This scaling is robust against some variations of the pump and of the lattice spacing.Furthermore, we compare the obtained scaling function with the one predicted for KPZ by functional Renormalization Group approaches, finding an excellent agreement.Finally, we provide compelling evidence for the 2D KPZ universality via the analysis of the fluctuations of the condensate phase.We compute the probability distribution function of the phase fluctuations, and show that it accurately matches the one obtained from large-scale numerical simulations of 2D KPZ systems [38,39,56].
The paper is organized as follows.The model studied, which is the generalized stochastic Gross-Pitaevskii equation for the dynamics of the condensate wavefunction, is introduced in Sec.II.We recall the mapping from this equation to the KPZ equation in Sec.III, and review the universal properties of the KPZ equation in 2D in Sec.IV.We describe the numerical simulations we performed, and we develop an analysis of the purely spatial and spatiotemporal vortices in Sec.V. We present our results for the spatiotemporal correlation function in Sec.VI, and our results for the distribution of the phase fluctuations in Sec.VII.

II. MODEL FOR THE POLARITON CONDENSATE DYNAMICS
The condensate wavefunction is described in 2D by a classical field ψ ≡ ψ(t, r) at time t and position r = (x, y).The condensate is populated through an incoherent excitonic reservoir of density n r (t, r) which is continuously driven by an external pump at a rate P .The reservoir density decays into the condensate via stimulated relaxation at rate R|ψ| 2 and into other modes at rate γ r .Under the assumption that the timescales for the reservoir and condensate dynamics are very different, the rate equation for the reservoir can be integrated out, which leads to the generalized stochastic Gross-Pitaevskii equation for the condensate [57,58] where m is the polariton effective mass, V is the lattice potential, g the polariton-polariton interaction strength, g r the repulsive interaction with excitons from the reservoir, n s = γ r /R, p = P R/(γ l,0 γ r ) = P/P th , with P th the threshold pump rate for condensation.We consider a kdependent loss rate of polaritons where k is the polariton in-plane momentum [18,59].Such a momentum dependence is equivalently introduced via a frequency-selective pump in the literature [52,60,61].The noise ξ is a complex variable with zero mean value, covariance ξ(t, r)ξ * (t , r ) = 2σδ(r − r )δ(t − t ), and an intensity set by the drive and losses as σ = γ l,0 (p + 1)/4.
As was shown in several studies, e.g.[14,15], the dynamics of the phase of a spatially homogeneous drivendissipative condensate can be mapped onto the KPZ equation.To show this, one expresses the condensate field in the density-phase representation ψ(t, r) = n(t, r)e iθ(t,r) .Under the assumption that the density fluctuations δn(t, r) = n(t, r) − n are sufficiently smooth ∇δn 0 and stationary ∂ t δn 0, one obtains for the phase dynamics the KPZ equation where the noise η is real and has zero mean and covariance η(t, r)η(t , r ) = 2δ(r − r )δ(t − t ).The KPZ parameters are determined by the microscopic parameters of Eq. ( 1) as ν = γ l,2 /2 + a /(2m), λ = − /m + aγ l,2 , D = σ(1 + a 2 )/(2n s (p − 1)), where a = 2pgn s 1 − 2P g r /(p 2 gγ r n s ) / γ l,0 .Several remarks are required.First, the mapping is readily extended to the case of an external confinement see eg. [16].Then, let us notice that for classical interfaces, the height field h(r, t) describes a D-dimensional manifold embedded in a d = D + 1 space.The mapping hence relates the phase θ(t, r) of a condensate in D = 2 spatial dimension to an interface of internal dimension two in a d = 2 + 1 space.It holds in any dimension provided the previous assumptions are fulfilled.Second, another fundamental difference lies in the compactness of the phase, which is periodically defined θ ∈ (−π, π], whereas the height of an interface is an unbounded field h ∈ (−∞, ∞).This implies that the phase field can feature topological defects, such as phase jumps or vortices as already mentioned.
While a compact version of the KPZ equation turns out to be relevant in some systems such as driven vortex lattices in disordered superconductors [62] or polar active smectic phases [63], we show that the behavior of the phase in the regime studied in this work belongs to the non-compact KPZ universality class, as it involves only a few sporadic vortices and the condensate phase can be uniquely unwound.As a third remark, let us mention that several works introduce some anisotropy in the polariton effective mass (m , m ⊥ ), which leads through a similar mapping to the anisotropic extension of the KPZ equation [13].However, as long as m and m ⊥ have the same sign, the anisotropy is irrelevant in the Renormalization Group sense and the universal properties of the anisotropic KPZ dynamics are controlled by the isotropic fixed point.When they have opposite sign, the non-linearity g KPZ = λ 2 D/ν 3 becomes irrelevant and the effective properties of the system are governed by the anisotropic Edwards-Wilkinson (Gaussian) fixed point [21,64,65].Since we are interested in the nonlinear KPZ regime, we do not include any anisotropy, and consider an isotropic polariton dispersion relation.
On the analytical side, the KPZ phase is described by a genuinely strong-coupling fixed-point in D ≥ 2. Indeed, it was shown in [69] that the perturbative RG flow equation for the KPZ coupling g KPZ can be resummed exactly, yielding an expression valid to all orders in perturbation theory in the vicinity of 2D.However, this expression fails to capture the strong-coupling KPZ fixed-point, yielding instead a run-away flow to infinity.Hence, let us emphasize that one cannot extract any quantitative information about the KPZ fixed-point in 2D or higher dimensions from the perturbative flow equation.
A non-perturbative RG calculation was devised using functional renormalisation group in Refs.[19,20,70], which allows accessing the strong-coupling KPZ fixed-point in all dimensions.It shows that the KPZ fixed-point is not connected to the Gaussian Edwards-Wilkinson fixed-point in any dimension, which explains the failure of perturbation theory, even resummed at all orders.The KPZ universal scaling function F in 2D was calculated in [20] and will serve as our theoretical reference.It is defined as where C is the two-point connected correlation function and C 0 and y 0 are non-universal normalisation constants.
The scaling function has the following asymptotics, where F 0 and F ∞ are constants given in Appendix C.

V. NUMERICAL SIMULATIONS AND ANALYSIS OF VORTICES
In the deep lattice regime, we solve a discretized version of Eq. (1) (see Appendix A for derivation and details) We numerically integrate the gGPE (6) using a split-step procedure for each noise realization (see Appendix A for details).We use a 2D N × N lattice grid with spatial spacing dx = 2.83µm, with N = 64 for the analysis of the correlation functions, and N = 32, 64, 128 for the study of the fluctuations distribution.In [18], where KPZ was observed in a 1D lattice, the full dispersion of the system was approximated in the simulations by a parabolic dispersion, which provided a very accurate description.
In 2D, we also approximate the dispersion of the lattice model (e.g. a cosine dispersion in the tight-binding regime) with a parabolic one.Indeed, the two dispersions are equivalent in the vicinity of k = 0 and differ only at large k.Thus, we expect the difference between the two to occur mainly at small spatial scales, which lie below the ones relevant for KPZ dynamics, and which are not of direct interest in this study.We choose the following parameters: effective mass m * = 8 × 10 −5 m e , where m e is the mass of the electron, corresponding to ReK/ = 0.72ps −1 ; γ l,0 = 0.31ps −1 , γ r = γ l,0 /10, γ l,2 = 0.1 µm 2 ps −1 , yielding ImK/ = 0.05ps −1 .We note that the value of γ l,0 corresponds to a polariton lifetime τ = 3.2ps, which is much too fast for significant equilibration to take place.Furthermore, we consider a non-interacting system with g = 0, since this turned out to be a good description of the 1D experiment.We note that the presence of polariton-polariton (g) and polariton-reservoir (g r ) interactions are not fundamental for the existence of a KPZ regime, since the mapping still holds for g = g r = 0.In fact, the study reported in Ref. [54] suggests that the development of the KPZ regime is favored by weak interactions.To further reduce the parameter space to explore, we also set g r = 0, although we stress that in order to faithfully model the experimental conditions one should take a non-zero value for this term.The saturation density is fixed at N s = 120.The numerical integration is performed with a time step of dt = 0.16ps to study the scaling, and down to dt = 0.0016ps to study the statistics of the fluctuations.This finer time discretization ensures that the phase can be uniquely unwrapped in time θ(t, r) ∈ (−π, π] → (−∞, ∞).The phase is unwrapped at a fixed space-point by constraining the phase difference between two consecutive times to be less than 2πζ.
In practice, we choose ζ not exactly one but ζ 1 due to the final time resolution.This procedure is robust as long as ζ is typically greater than 0.5.
In order to tentatively locate the KPZ regime, we first investigate the density of topological defects, which are either purely spatial vortices in the (x, y) plane, or spacetime vortices which have a non-vanishing projection of their vorticity in the (t, x) and/or (t, y) planes.The latter were shown to play an important role in 1D [17,18].Maintaining the other parameters fixed, we vary the reduced pump p from low p 1 (close to threshold) to moderately high p = 2.For each value of p, we determine the density of purely spatial and space-time vortices in the stationary state, which is shown in Fig. 1, with typical phase configurations in the (x, y) plane (similar maps are obtained in the (t, x) and (t, y) planes).The numerical procedure to detect the vortices and compute their density is detailed in Appendix A. We find that the number of vortices drastically decreases when increasing the pump power, and for p 1.6, very few spatial vortices are present, in agreement with Ref. [52], and also very few space-time vortices.Both appear only by pairs of closely located vortex and anti-vortex.We hence focus in the following on the highest value of the pump p = 2, which appears as more favorable for observing KPZ dynamics.We show in Appendix B that the KPZ regime is robust, to a certain extent, against decreasing the pump power.
The presence of topological defects turns out to be also sensitive to the discreteness of space, and we find that the presence of a lattice (of spacing dx) is favorable to observe the KPZ regime since the density of such defects is suppressed when increasing the grid spacing dx, as shown in Fig. 2. For dx 1µm, the emergence of space-time vortices is rare and very short-lived.In the following, we choose dx = 2.83µm, which is close to the spacing in experimental micro-pillar lattices [18].For this value, the fraction of both spatial and space-time vortices is below 10 −3 and does not spoil the KPZ universal properties.Moreover, the core size of the spatial vortices was observed to be very narrow, of the order of the lattice spacing.
modulus ∆r = |∆r|, therefore we also average the data over circular shells of radius ∆r around r = 0.Under the assumption that the density-phase and density-density correlations are both negligible, the expression of g 1,ψ becomes g 1,ψ (∆t, ∆r) e i∆θ , where ∆θ ≡ θ(t 0 + ∆t, r 0 + ∆r) − θ(t 0 , r 0 ).Upon performing a cumulant expansion of g 1,ψ to lowest order, one deduces that it is related to the connected correlation function C of the phase as where C is the analogue quantity as Eq. ( 4) for the in-terface.Note that g 1,ψ constitutes a reliable observable to access the scaling behavior of the phase only provided the three previous assumptions are verified.This is not guaranteed a priori and has to be assessed.It was shown to be the case in 1D in the experimental conditions where the KPZ regime was evidenced [18].They are also satisfied in our study, as shown in Appendix B.
We first study the equal-time and equal-space correlation functions.If the phase follows the KPZ dynamics, then according to Eqs. (3, 5) with χ 0.39 and β 0.24.We observe the expected scaling behavior both in space and time, as shown in Fig. 3.For temporal correlations, the KPZ scaling extends for time differences spanning over a decade, from 3 × 10 2 − 1.5 × 10 4 ps, whereas for spatial correlations, the KPZ scaling is observed over a range from 10 − 40µm.The limited space range is related to the relatively small size of our system, with a condensate of approximately 90µm radius.Note that the polariton interaction with the reservoir g r was here set to zero, but it is likely to play a role on the typical space and time scales where KPZ dynamics dominates, as observed in 1D [18].
In order to provide a quantitative estimate of the critical exponents χ and β, we compute the local exponents, given by the following logarithmic derivatives The result is shown in the insets of Fig. 3.We fit the obtained local exponents by a constant in the appropriate spatial and temporal windows, which yields the values of the universal exponents β EP 0.22 ± 0.06 and χ EP 0.36 ± 0.04.These values are in good agreement with the results from numerical simulations for the 2D KPZ universality class [35].
We now study the scaling properties of the spatiotemporal correlations over the whole relevant domain (∆t, ∆r).For this, we first select the data lying within the scaling regime, indicated in the inset of Fig. 4. We then construct the scaling function F (y) defined in Eq. ( 3).The normalizations C 0 , y 0 are determined from fitting the equal-time and equal-space correlation functions with the power-laws Eqs.(9a, 9b) in the appropriate ranges.Our results are shown in Fig. 4, together with the theoretical scaling function calculated from functional Renormalization Group.We observe a collapse of the data onto a single curve, which matches the theoretical one with excellent accuracy.This result shows that the phase of the two-dimensional polariton condensate follows a KPZ effective dynamics over extended length and time scales.

VII. STATISTICAL PROPERTIES OF PHASE FLUCTUATIONS
We next compute the probability distribution of the fluctuations of the unwrapped phase θ ≡ θ(t, r 0 ) at a fixed space point r 0 .In order to discuss the 2D results, we first briefly summarize what is known for a 1D classical interface described by a height field h.In this case, the distribution of the reduced height fluctuations defined as h = (h − v ∞ t)/(Γt) 1/3 , with v ∞ and Γ non-universal parameters, is known exactly.In fact, the distribution of h was shown to be universal but sensitive to the global geometry of the growth, which led to the definition of three main universality sub-classes.For flat, curved, or stationary geometries, the probability distribution is given by Tracy-Widom GOE (Gaussian Orthogonal Ensemble), Tracy-Widom GUE (Gaussian Unitary Ensemble) or Baik-Rains distributions, respectively [71].They all share the common feature of being non-Gaussian, with finite skewness and excess kurtosis.Strikingly, the three universality sub-classes can also be realized in the polariton condensate in 1D, where the role of global geometry can be emulated by appropriate external potentials [16].
In 2D, like in 1D, the height field exhibits a linear growth in time with average velocity v ∞ , over which fluctuations develop with a t β scaling.All 2D KPZ numerical works [38,39] find that the distribution of the reduced fluctuations h are non-Gaussian, and suggest, like in 1D, the existence of different geometry-dependent universality subclasses.However, the corresponding distribution functions have no known analytical form like Tracy-Widom or Baik-Rains distributions in 1D.To in- (fuchsia, light green, dark gray symbols, respectively).Each curve comprises of approximately 8 × 10 6 data for each system size.The centered and unit-variance Gaussian distribution (denoted N (0, 1)) is also shown for comparison, together with the distribution found from KPZ large-scale simulations (denoted KPZ num.) [39].
vestigate the statistical properties of the phase fluctuations, we record 5120 independent realisations of the time evolution of ψ, from which we extract the phase and unwrap it in time at a fixed space point.The presence of a few space-time vortex-anti-vortex pairs with random core positions induces phase slips at the position of interest, which are non-instantaneous phase jumps close to (but lower than) 2π.
As already evidenced in 1D [18], these phase slips result in the existence of replicas of the distribution of the fluctuations separated by approximately 2π.As detailed in Appendix D, one can select the main central distribution at each time instants in the KPZ window and sum them after appropriate re-scaling.The obtained distribution is shown in Fig. 5.It is markedly distinct from a Gaussian distribution of same mean and variance, and matches the distribution [72] computed in 2D KPZ large-scale numerical simulations for flat interfaces [38,39], which is the expected universality sub-class for polaritons in the absence of confinement potential [16].This result for the fluctuations distribution enforces the evidence for KPZ universality in 2D polariton condensates.

VIII. CONCLUSIONS AND PERSPECTIVES
In this work, we have shown, using numerical simulations of the discrete stochastic generalised Gross-Pitaevskii equation, that a KPZ regime can be achieved in discrete 2D driven-dissipative exciton polariton condensates under incoherent pumping.
We have obtained the condensate spatiotemporal first-order correlation function, and shown that, for the parameters studied, it exhibits stretched exponential behavior, both in space and time, with critical exponents characteristic of the 2D KPZ universality class.This scaling persists in a finite region of pump strength and grid spacing.Moreover, the associated scaling function accurately matches with the theoretical KPZ scaling function given by functional Renormalization Group methods.We have also obtained the distribution of the phase fluctuations, which is highly non-Gaussian and very close to the distribution computed in numerical simulations of KPZ interfaces.This is a compelling evidence that the phase fluctuations behave as the KPZ stochastic process.
Our findings open promising perspectives.On the theoretical side, it would be desirable to obtain a complete description of the phase diagram of a polariton condensate in 2D, with a clear understanding of the interplay between the various possible regimes, as e.g.nonequilibrium BKT, and KPZ.This remains a challenge as the parameter space to explore is very large.On the experimental side, the demonstration of the KPZ scaling in a 1D polariton condensate was recently realized in a lattice, where real space was discretized [18].Since our simulations predict that this is the optimal situation to prevent topological defect proliferation in 2D, similar techniques could be used in 2D.The evidence of a KPZ regime in a 2D polariton condensate would constitute a major breakthrough, especially since a convincing experimental realization of KPZ universality class in 2D is still missing and actively sought for.
Appendix A: Discrete gGPE, numerical integration, detection of vortices Our work is based on the numerical solution of the generalized Gross-Pitaevskii equation ( 1) by discretizing it on a grid of spacing dx.This may describe both a continuous system, by taking the limit of lattice spacing tending to zero, as well as a specific lattice geometry, which could be realized in the polariton experiment.We provide in this Appendix some further information concerning the mapping to the lattice model, the numerical integration and the method for identifying the vortices.

Derivation of the lattice model
In order to describe the case of polaritons on a lattice, we start from Eq. ( 1) where V (x) is the periodic lattice potential.In the deep lattice regime, the condensate wavefunction can be expanded on the Wannier function basis ψ(t, x) = j w j (x)ψ j (t), with w j (x) = w(x − x j ) where w(x) is localized around x = 0.The Wannier functions are a set of real, localized, orthonormal wavefunctions centered on the lattice sites x j , where j = {j x , j y } identifies each lattice site of the 2D grid.
We shall assume in the following for simplicity that the lattice potential is separable in the x, y directions, i.e. can be written as V (x) = V (x) + V (y), and we denote by dx the lattice spacing.Simple examples of V (x) are the sinusoidal case V (x) = V 0 sin 2 (k 0 x), with k 0 = π/dx or the Kronig-Penney potential V (x) = V 0 n Θ(x − (ndx + b))Θ((n + 1)dx − x), with b being a parameter associated to the width of the potential barrier and Θ(x) the Heaviside step function.If the potential is separable, one also has w(x) = w 1d (x)w 1d (y) with w 1d (x) the Wannier function for the one-dimensional lattice problem.We will further assume that |ψ| 2 /n s 1 such that the terms involving the pump can be Taylor expanded to first order.This hypothesis may be easily released by treating in an analogous fashion the higherorder terms of the Taylor expansion.
Under the above assumptions, by projecting Eq. (A1) onto the Wannier function basis and by neglecting overlaps among Wannier functions on different sites in the non-linear terms, we obtain the following discrete non-linear Schrödinger equation i ∂ t ψ j = − K(ψ jx+1,jy + ψ jx−1,jy + ψ jx,jy+1 + ψ jx,jy−1 ) where the parameters entering in Eq. (A2) are In the numerical simulations, in order to make link with the continuous limit expressions, we have resummed the term (p(1 − |ψ j | 2 /N s ) onto p/(1 + |ψ j | 2 /N s ).Also, for convenience, we have chosen the zero of the energy such that W 0 = 4K.
The continuous limit of Eq. (A2) can be formally defined by the limit ψ(x j ) = lim dx→0 ψ j /dx.This allows one to relate the parameters of the lattice model to the ones of the effective continuous model, for example Re(K) = 2 /(2m * dx 2 ), with m * an effective mass.However, this limiting procedure assumes that the parameters entering the discrete GPE are independent from the lattice spacing, which is -strictly speaking-not the case since the lattice potential enters in the definition of K.This implies that the effective mass will be slightly modified when changing dx.

Details of the numerical integration
The integration of the discretized Eq. ( 1) is performed using the split-step Fourier method.This scheme consists of treating the kinetic and potential/interaction parts of the equation separately when propagating in time.More specifically, after suitable discretization in space, the wavefunction after a single time step dt is constructed by the following three building blocks: i) propagation by dt in real space using the potential and interaction parts, ii) propagation of dt in momentum space by the kinetic part, and iii) addition of the stochastic contribution (Wigner noise).This procedure is repeated in a loop until the desired final time is reached.We note that, switching back and forth from real to momentum space is done particularly efficiently via the Fast Fourier Transform procedure.The source code of our solver is available at this repository.and the one from the phase correlations g 1,θ (∆t, ∆r) = e i(θ(t0+∆t,r0+∆r)−θ(t0,r0)) .(B2) The results for p = 2 and p = 1.6 are displayed in Fig. 7.We observe that g 1,n is very close to unity and constant over the whole KPZ spatial and temporal ranges for both pump powers.Thus, the role of density-density correlations is indeed negligible.The density-phase time correlations are almost vanishing, since g 1,ψ (∆t, 0) and g 1,θ (∆t, 0) almost perfectly coincide over the whole KPZ time range, as shown in panel (a) of Fig. 7.The densityphase space correlations are more important, and their effect increases while decreasing the pump.They are shown in panel (b) of Fig. 7.For p = 2, the KPZ stretched exponential well models the behavior of both g 1,ψ and g 1,θ .We checked that the difference between the two curves is just a multiplicative factor.Hence, the density-phase correlations are almost constant in the spatial KPZ window, and the scaling of g 1,ψ follows the scaling of the phase-phase correlations, as expected in a KPZ regime.For p = 1.6, the KPZ stretched exponential no longer provides a good fit for the behavior of g 1,ψ (as already visible in Fig. 7) and neither for g 1,θ .Moreover the density-phase correlations start varying within the KPZ window, the difference between g 1,ψ and g 1,θ is no longer a mere multiplicative factor.Thus, one concludes that for p = 2, the behavior of the correlations g 1,ψ indeed reflects the KPZ scaling of the phase itself, both in space and time, whereas for p = 1.6, the density-phase space correlations also contribute to the behavior of the g 1,ψ , and the KPZ regime is hindered in space, although it persists in time.
at zero, by subtracting the first moment of the selected data, δθ ti = θ ti − θ ti , and normalized to unit variance.We then sum the fluctuations of the first peaks, properly normalized, for all time instants in the appropriate KPZ time window.Let us mention that we have also performed a detailed analysis of the skewness and excess kurtosis of the phase fluctuations.These quantities are insensitive to normal-isation issues, whereas the data analysis of the distribution can suffer from unavoidable errors due to the empirical choice of and P min , and the small number of discrete bins for given t i .This analysis is reported in [74], and confirms that the skewness and excess kurtosis computed from our data is in good agreement with the universal values obtained from KPZ numerical simulations [38,39].

2 FIG. 1 .
FIG.1.Density of purely spatial and space-time vortices as a function of the pump power p, with typical phase configurations in the (x, y) plane (vortex -green, anti-vortex -red) shown for p = 1.2 (left) and p = 1.8 (right) for a grid spacing of dx = 2.83µm.

FIG. 2 .
FIG.2.Density of purely spatial and space-time vortices as a function of the lattice spacing, with typical phase configurations in the (t, y) plane (vortex -green, anti-vortexred) shown for a grid spacing of dx = 0.71µm (left) and dx = 1.41µm (right), for a fixed pump power p = 2.

FIG. 3 .
FIG. 3. (a)Equal-space and (b) equal-time correlation functions (blue plain lines) evaluated from g 1,ψ as explained in the text, compared with the KPZ scaling laws (black dashed lines).The numerical data follow the KPZ scaling for an extended range of time and space separations (indicated by the gray shade).The local exponents, which are computed from Eqs. (10a, 10b) are shown in the inset in both cases, together with their fits by a constant in the gray region (solid green lines, denoted 2βEP and 2χEP ) and the values from KPZ numerical simulations (black dashed line, denoted 2β and 2χ).Note that a running average has been used in order to smoothen the equal-space correlation function (i) before computing its derivative.

FIG. 4 .
FIG.4.Universal scaling function F (y) from our numerical data (blue dots) and from the functional Renormalization Group calculation (black solid line).Inset: Space-time map of g 1,ψ , with the red contour delimiting the scaling region where data are selected to construct F (y).

FIG. 5 .
FIG.5.Histogram of the centered and rescaled to unitvariance phase fluctuations δθ in the appropriate KPZ time window for different system sizes L = 90.5µm,181µm, 362µm (fuchsia, light green, dark gray symbols, respectively).Each curve comprises of approximately 8 × 10 6 data for each system size.The centered and unit-variance Gaussian distribution (denoted N (0, 1)) is also shown for comparison, together with the distribution found from KPZ large-scale simulations (denoted KPZ num.)[39].
FIG.6.(a) Equal-space and (b) equal-time correlation functions g 1,ψ , with darker colours corresponding to larger values of the pump parameter, p = 1.6, 1.7, 1.8, 2, 2.5.The correlations are compensated by the appropriate KPZ power laws, such that the corresponding data appear as plateaus in the KPZ spatial and temporal ranges, which are highlighted in gray.

FIG. 7 .
FIG. 7. (a) Equal-space and (b) equal-time correlation functions g 1,ψ (red) and g 1,θ (blue), for p = 1.6 (lighter shade) and p = 2 (darker shade).The dashed lines indicate the corresponding stretched exponentials with KPZ exponents.In panel (b), the lighter shades are shifted upwards after multiplying by a factor of 1.25 in order to avoid overlapping with the darker shades.Inset: Density-density correlation g1,n for these values of the pump.