{'fraction_non_alphanumeric': 0.045661363709703445, 'fraction_numerical': 0.01445674271515699, 'mean_word_length': 4.344946533287341, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Space-time vortex driven crossover and vortex turbulence phase transition in one-dimensional driven open condensates

We find a first order transition driven by the strength of non-equilibrium conditions of one-dimensional driven open condensates. Associated with this transition is a new stable non-equilibrium phase, space-time vortex turbulence, whose vortex density and quasiparticle distribution show strongly non-thermal behavior. Below the transition, we identify a new time scale associated with noise activated unbound space-time vortices, beyond which the temporal coherence function changes from a Kardar-Parisi-Zhang type subexponential to a disordered exponential decay. Experimental realization of the non-equilibrium vortex turbulent phase is facilitated in driven open condensates with a large diffusion rate.

The last decade has witnessed fast experimental development in realizing driven open quantum systems with many degrees of freedom. Examples include excitonpolaritons in semiconductor heterostructures [1][2][3], ultracold atoms [4][5][6], trapped ions [7,8], and microcavity arrays [9,10]. The common characteristic is explicit breaking of detailed balance on a microscopic level by the presence of both coherent and driven-dissipative dynamics, placing these systems far from thermal equilibrium. This makes them promising laboratories for studying non-equilibrium statistical mechanics, according to which one expects non-equilibrium features to persist to the macroscopic level of observation.
A case in point are exciton-polariton systems, which can be engineered in one-and two-dimensional geometries [1][2][3]11]. Due to effectively incoherent pumping, these systems possess a phase rotation symmetry and can thus show Bose condensation phenomena. While their dynamics is described microscopically in terms of a stochastic complex Ginzburg-Landau equation (SC-GLE) [12], recently it was noticed that at low frequencies it maps to the Kardar-Parisi-Zhang (KPZ) equation [13,14]. Traditionally, the latter equation describes, e.g., the roughening of surfaces, with the dynamical surface height being unconstrained [15]. In contrast, the dynamical variable in the context of driven open quantum systems is the condensate phase, which is compact. The compact KPZ (cKPZ) equation gives rise to a novel scenario for non-equilibrium statistical mechanics realized in concrete experimental platforms, and raises the fundamental question of the physical consequences of compactness.
In this work, we address this question for onedimensional driven open quantum systems. To this end, we establish the complete phase diagram of 1D driven open condensates (DOC) via numerical simulations in combination with analytic analysis, and explain its basic structure on the basis of the 1D cKPZ equation (cf. Fig. 1). It can be traced back to the behavior of dynamical topological defects, namely space-time vortices ( Fig. 1(b)). More specifically, we find: (i) The emergence of a new time scale t v in the long time behavior of the temporal coherence function. At weak noise level σ, it is exponentially large, t v ∝ e B/σ (B a non-universal positive constant), reflecting an exponentially suppressed but finite space-time vortex density ( Fig. 1(c)). The subsequent asymptotic regime t ≫ t v is characterized by a disordered, exponentially decaying first order temporal coherence function. The exponential dependence on the inverse noise level corroborates previous results on the observability of KPZ physics in 1D, where the crossover scale from stretched exponential equilibrium diffusive to non-equilibrium KPZ scaling behaves algebraically as t * ∝ σ −2 [16] at weak noise level, ensuring generically t * ≪ t v [17]. (ii) We identify a new intrinsic non-equilibrium phase of the cKPZ equation, with the following key signatures: (a) A first order transition at low noise level from a regime of exponentially low to high space-time vortex density (cf. Fig. 1(d)). (b) Strongly non-thermal behavior of the momentum distribution of quasiparticles n q ∝ q −γ at large momentum q (cf. Figs. 3(a) and 3(b)), with γ significantly deviating from the value for thermal behavior (γ ∼ 2). Such a strong scaling behavior is reminiscent of turbulence, which usually is a transient phenomenon that occurs in specifically initialized systems undergoing purely Hamiltonian dynamics without external drive [20]. In contrast, it occurs in stationary state in our case.
Microscopic model.-We describe the dynamics in terms of the SCGLE with complex Gaussian white noise of zero mean, as appropriate for experiments with excitonpolariton systems [1][2][3]. It reads in 1D [12] (units = 1) . r c and u c are the chemical potential and the elastic collision strength, respectively. r d = γ p − γ l is the difference between the single particle loss γ l and incoherent pump γ p . For the existence of a condensate in the mean field steady state, r d > 0. u d is the positive two-particle loss rate; K c = 1/(2m LP ) with m LP the effective polariton mass and K d a diffusion constant. We obtained most of the numerical results presented in this work by solving Eq. (1) using the same approach as in Ref. [21] (see also [22] for technical details) and we set r d = K d = 1; hence, t and x are measured in units of r −1 d and √ K d . If not specified otherwise, we used 10 3 stochastic trajectories to perform ensemble averages. Low frequency effective description.-In the absence of phase defects, the low frequency dynamics of the system is effectively described by the KPZ equation [15] for the phase of the condensate field, [13]. Here, D describes phase diffusion, and ξ is a Gaussian white noise of strength 2σ KPZ . The non-linearity λ is a direct, single-parameter measure for the deviation from equilibrium conditions, with λ = 0 in the presence of detailed balance [13]. In order to properly account for the compactness of the phase and to allow for a description of phase defects, we work with a lattice regularized version of the KPZ equation (cKPZ), which can be straightforwardly derived from a spatially discretized SCGLE on a 1D lattice with spacing ∆ x [22]. It reads is reproduced upon assuming that phase fluctuations are small, and taking the continuum limit. The crucial difference between the non-compact continuum KPZ and the compact KPZ equation is revealed by the number of independent scales in the problem, which originates from the compactness of the phase. Indeed, by rescaling t, ξ, and θ in the continuum case, there is only one tuning parameter given by g ≡ λ(σ KPZ /2D 3 ) 1/2 . In contrast, for the cKPZ equation, we can rescale t andξ i but not the phase field θ i due to its compactness, resulting in two independent tuning parameters. Rescaling amounts to replacingD → 1,λ →λ = λ/D,σ →σ = ∆ x σ KPZ /D in Eq. (2). The rescaled equations withλ and −λ are equivalent, therefore we further redefineλ ≡ |λ/D|. The fact thatλ andσ are two independent tuning parameters suggests that there must exist new physics associated with changing each of them, beyond the physics of the KPZ equation, where changing the nonequilibrium strengthλ and noise strengthσ are equivalent to changing the single parameter g. Indeed, as we shall see in the following, the noise strength is associated with a new time scale in the 1D DOC giving rise to an additional scaling regime in dynamical correlation functions. Moreover, the nonequilibrium strengthλ can drive the system to a new non-equilibrium vortex turbulent phase via a first order transition at low noise level.
Space-time vortex driven crossover forλ/λ * < 1.-The prime observable that distinguishes a 1D DOC from its equilibrium counterpart [21,26] is the autocorrelation function C ψ t (x; t 1 , t 2 ) ≡ ψ * (x, t 1 )ψ(x, t 2 ) . Its long time behavior is determined by fluctuations of the phase, C ψ t (x; t 1 , t 2 ) ∝ e −c∆ θ (t1−t2) , where c is a non-universal constant, and ∆ θ ( [27], with L being the spatial size of the system. Previous numerical studies [21,26] in the regime of weak noise, defined by the absence of space-time vortices (cf. Fig. 1(b)) in the spatio-temporal extent of the numerical experiments. Increasing the noise level leads to proliferation of space-time vortices, which in turn strongly affect the temporal coherence of DOCs as we describe in the following.
Figure 2(a) shows ∆ θ (t 1 − t 2 ) for moderate noise strengths. For the lowest value σ = 8 −1 , KPZ scaling can be observed in a wide range 5 × 10 3 |t 1 − t 2 | 5 × 10 4 , after which ∆ θ (t 1 − t 2 ) grows linearly with time. Accordingly, the autocorrelation function C ψ t (x; t 1 , t 2 ) exhibits a crossover from stretched-exponential to simple exponential decay at long times [22]. At weak noise, the crossover time t c , defined as the point where the gradient of ∆ θ (t 1 − t 2 ) on a double logarithmic scale exceeds 0.9, increases exponentially with the inverse noise level, i.e., t c ∝ e C/σ , where C is a positive constant. This is shown in Fig. 2(b). The fast growth of phase fluctuations and the associated decoherence of DOCs for strong noise and at long times is due to space-time vortices. Numerical evidence for this connection is presented in Fig. 1(c), where a linear fit on the semi-logarithmic scale clearly demonstrates that the space-time vortex density behaves as P v ∝ e −A/σ (see [22] for details on the numerical cal- culation of P v ) in the range 6 σ −1 13, while saturation sets in at even higher noise strengths. This exponential determining P v is the same as the one appearing in t c and can be interpreted as the statistical weight of space-time vortex configurations in an equivalent 1+1D static equilibrium description of the problem. Formally, this 1+1D static equilibrium description can be established by rewriting the cKPZ equation as a functional integral over space-time configurations of the phase field as Z =´D[θ]e −H[θ]/T [22,29], with the Boltzmann factor e −H[θ]/T determining the probability with which a particular solution θ(x, t) is encountered in a numerical integration of Eq. (2) from time t 0 to t 1 in a system of size L. The temperature is T = 4σ KPZ , and the effective Hamiltonian (going back to a spatial continuum notation) is . Regarding z = t as a spatial coordinate, this is the Hamiltonian of a 2D smectic A liquid crystal [29]. Crucially, one can show that the energy [30], indicating that at any finite noise strength there is a non-zero density of noise activated vortices P v ∝ e −Ev/T . Deviations from this exponential dependence in Fig. 1(c) for strong noise are due to the interaction energy of vortices, which becomes important at higher densities.
The exponentially small but finite space-time vortex density entails a time scale t v [32] characterizing the average temporal distance between vortices, t v ∝ P . Then, using a simple phase random walk model, one can show that space-time vortices occurring at multiples of t v lead to linear growth of ∆ θ (t 1 −t 2 ), hence disordered behavior of the phase correlations, beyond a time scale t c = O(t v ) [22].
First order transition and vortex turbulence forλ/λ * > 1.-We now investigate the strongly non-equilibrium regime at largeλ, where most of the results presented in the following are obtained from direct simulations of the SCGLE. Figure 1(d) shows the dependence of the vortex density P v onλ at weak noise (σ = 10 −2 ). At smallλ , P v is exponentially small in line with the discussion above. However, when tuning above a critical strength λ * ≃ 3.23, the vortex density undergoes a sudden jump by around 10 orders of magnitude, indicating a sharp first order transition at low noise level [35].
Clear signatures of the transition can also be seen in the momentum distribution function n q ≡ ψ * (q)ψ(q) , which is accessible in experiments with exciton-polaritons. At large q, it behaves as n q ∝ q −γ (cf. Figs. 3(a) and 3(b)). Across the transition the value of γ undergoes a jump from γ ≃ 2, which is characteristic of noise activated vortices [20], to γ ≃ 5 (with a weak dependence on parameters). Such strong scaling behavior is reminiscent of turbulence and we thus refer to this phase as vortex turbulence (VT).
The physical origin of the VT phase and the associated first order transition can be traced back to a dynamical instability triggered byλ: consider the dynamical equations for the phase differences between nearest neighboring sites, The first term, originating from the diffusion term in the cKPZ, is a restoring force which attenuates phase differences, while the second term amplifies them. This causes the dynamical instability for largẽ λ and induces vortices without resorting to noise induced excitation. The dynamical instability is thus triggered on short scales. Moreover, this mechanism can also occur in higher dimensions. Taking into account the exponential suppression of noise activated vortices at low noise and the fact that the vortex generation due to the dynamical instability is insensitive to weak noise, this rationalizes the existence of a first order transition tuned byλ. Increasing the noise strength, Figs. 3(c) and 3(d) show that λ * does not depend on σ within the numerical accuracy of our simulations (λ * = 3.23 with an error bar of ±0.07). The transition line terminates at σ * ≃ 0.014 at an apparent second order critical point. It is characterized by a vanishing jump ∆P v , whose derivative with respect to σ diverges according to ∝ (σ − σ * ) −κ with κ ≃ 0.63.
The discussion of the VT transition reveals it to be rooted in the non-equilibrium nature of the cKPZ equation [36]. Numerically solving this equation at low noise level, we obtainλ * ≃ 20. This can be compared to numerical simulations of the SCGLE at weak noise, which yield the significantly reduced valueλ * ≃ 3.23. The discrepancy is due to the mutual feedback between phase and density fluctuations present in SCGLE: the former can cause the latter via a phase-density coupling term proportional to the diffusion constant K d (cf. Eq. (1)), which in turn facilitates strong phase fluctuations, gives rise to vortex creation and consequently causes the critical valueλ * in the SCGLE to depend on K d .
Experimental observability.-In exciton-polaritons [11], typically the diffusion constant K d ≪ K c = 1/2m LP (cf. Eq. (1)). Then, in line with previous results [21,26], KPZ scaling should be observable due to the exponential suppression of unbound vortices at low noise and weak non-equilibrium strength. However, we note that a relatively largeλ (but below the VT transition) is generally favorable in order to overcome the influences from potential inhomogeneities and the finite system size. On the other hand, K c ≃ 0 in recently realized excitonpolariton condensates formed in the flat band of a 1D Lieb lattice of micropillar optical cavities [38]. In this case, the maximum achievable non-equilibrium strength isλ max ∼ 2K d /K c , indicating thatλ can be tuned large to make the VT phase accessible. As explained above, the distinct features of this phase are revealed in the momentum distribution and can thus be obtained by momentum resolved correlation measurements [12]. Alternatively, explicit engineering of a large diffusion constant could be achieved in 1D arrays of microwave resonators coupled to superconducting qubits [39]. More generally, DOC systems with relatively large K d /K c are favorable to observe the VT phase, since the associated mutual feedback between phase and density fluctuations makes it easier to trigger the dynamical instability.
Conclusions.-The non-equilibrium phase diagram of one-dimensional driven open condensates is crucially shaped by space-time vortices, as the relation to the compact KPZ equation reveals: at weak non-equilibrium strength, they govern the asymptotic behavior of the temporal correlation functions, however only beyond an exponentially large crossover time scale. This protects KPZ physics and suggests its observability in current excitonpolariton experiments. Moreover, these defects cause the existence of a new phase under strong non-equilibrium conditions, stationary vortex turbulence. We believe that our predictions will stimulate both further theoretical research on the cKPZ equation, especially in higher dimensions and in the context of stationary turbulence, as well as experimental efforts in searching for these two non-equilibrium phases. Other intriguing directions for future research are the study of nucleation dynamics and phase coexistence in the vicinity of the transition to VT, and the critical behavior of the second order end point using the techniques developed in [19].
Note added.-Upon completion of this manuscript, we became aware of the work by Lauter et al.
[40], reporting a related dynamical instability in arrays of coupled phase oscillators.
We thank I. Carusotto, B. Kim, E. Altman, D. Huse, J. Toner and S. Mathey for useful discussions, and the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for allocation of computer time. This work was supported by German Research Foundation (DFG) through ZUK 64, through the Institutional Strategy of the University of Cologne within the German Excellence Initiative (ZUK 81) and by the European Research Council via ERC Grant Agreement n. 647434 (DOQS). L. S. acknowledges funding through the ERC synergy grant UQUAM. imizes this Hamiltonian subject to the topological constraint that the circulation of the phase around the location of the vortex core is an integer multiple of 2π. In Ref.
[31], the energy Ev was shown to be finite for λ = 0. When λ = 0, it is not possible to find θv(x, t) minimizing the energy analytically, but an upper bound on Ev can be obtained by calculating H[θv] with the λ = 0 vortex configuration. A straightforward integration shows that this upper bound is finite. This is in sharp contrast to a superfluid in 2D in thermal equilibrium for which the energy of a single vortex diverges logarithmically with system size.
[32] The rate 1/tv at which phase slips occur and hence the vortex density Pv for a system in thermal equilibrium (described by Eq. (1) with rc = Kc = uc = 0) was determined in the context of narrow superconducting channels in Refs. [33]. Phase slips are fluctuations which reduce the current and thus leads to a finite resistance of the channel. More formally, during a phase slip the system changes its state from one local minimum of the free energy to a neighboring one that differs in current by 2π/L. At low noise (set by the temperature T ), the rate at which such fluctuations occur is proportional to exp(−δF/T ) where δF is the height of the free energy barrier separating the local minima. Out of thermal equilibrium, however, the notion of a free energy becomes void and the theory developed in Refs.
[33] cannot be applied directly. Then, the approach presented in the main text is needed to obtain the exponential factor. Here we provide some details on the derivation of the compact KPZ equation in one spatial dimension, and the equivalent equilibrium description in terms of an effective Hamiltonian.

Compact KPZ equation
Mathematically, the compactness of the phase field can be formulated as a "discrete" local gauge symmetry of the SCGLE in the amplitude-phase representation of the field ψ(x, t) ≡ ρ(x, t)e iθ(x,t) . Namely, if ρ(x, t)e iθ(x,t) is a solution of the SCGLE, then ρ(x, t)e i[θ(x,t)+2πn(x,t)] is also a solution, where n(x, t) is an integer valued function. This symmetry reflects the simple physical requirement that any choice of the space-time local definition θ(x, t) by shifting 2πn(x, t) must not change the physical results. Therefore, it is crucial to require the effective description to respect this symmetry.
The discrete nature of the gauge transformation indicates that both the temporal and spatial derivative of the phase field have to be taken with care when one tries to derive the dynamical equation for the phase field from the SCGLE. This is due to the fact that the existence of the derivatives for any physical solution of the phase field is not always guaranteed. To see this point, let us consider a solution of the SCGLE, ψ(x, t) = ρ(x, t)e iθ(x,t) with θ(x, t) having well defined ∂ t θ and ∂ x θ at any space-time point. For a generic gauge transformed solutionψ = ρ(x, t)e iθ(x,t) withθ(x, t) = θ(x, t) + 2πn(x, t), one immediately sees that ∂ tθ or ∂ xθ is not well defined if around a certain space-time point, say ( . Nevertheless, due to fact that the exponential function e iα is a smooth function with respect to α for any physical solution of the phase field, one can always choose a gauge n(x, t) in such a way that either the temporal or the spatial derivative of the transformed phase field exists everywhere in the space-time domain under consideration, but not both.
In the following, we choose the gauge in which the temporal derivative of the phase field exists everywhere. We choose the value of n(x, t + 0 + ) according to θ(x, t) in such a way that θ(x, t + 0 + ) − θ(x, t) is infinitesimal. Then, the spatial derivative of the phase field is not well-defined everywhere. At the price of introducing a short distance scale, we can circumvent this issue by first rewriting the SCGLE on a discrete 1D lattice with lattice spacing ∆ x , such that the second order spatial derivative is replaced by the a second order finite difference. The discretized SCGLE reads The rest of the derivation goes along the lines of the similar derivation presented in [1]. Plugging the amplitude-phase decomposition of ψ, i.e. ψ i (t) = (M 0 + χ i (t))e iθ i (t) with M 0 the homogeneous mean-field solution, into (S-1), we arrive at a coupled discretized equation of motion (EOM) of amplitude and phase fluctuations, We linearize in χ and eliminate it adiabatically, made possible since its evolution is fast compared to the gapless phase degree of freedom. From the EOM of χ i we get Plugging the above equation into the EOM of θ i , neglecting sub-leading and non-linear amplitude fluctuations ∂ t χ i , χ 2 i , χ 3 i also here, and keeping only the additive part of the noise, we get the compact KPZ equation Here, we construct the effective 2D static equilibrium model for the compact KPZ equation. Conceptually, the idea is to interpret the compact KPZ equation as a mapping from the set of stochastic variablesξ(t) = {ξ i (t)} with Gaussian path probability distribution to new stochastic variables θ(t) = {θ i (t)}. The distribution of the new variables can be obtained using the usual relation [2] where T ≡ 4σ KPZ is the effective temperature and H[θ] assumes the form In deriving P θ , we used that the Jacobian determinant associated with the change of variables fromξ to θ is equal to unity for a retarded regularization of the time derivative in the compact KPZ equation. In the continuum limit, H[θ] reduces to the form shown in the main text. Finally, calculating the normalization of P θ corresponds to calculating the partition function Z of the model defined by Eq. (S-10) at temperature T ,

CROSSOVER BEHAVIOR IN THE TEMPORAL CORRELATION FUNCTION
A crossover behavior similar to the one seen in the temporal phase fluctuation ∆ θ (t 1 − t 2 ) is also observed in the modulus of the condensate temporal auto correlation function In practice, by assuming spatial translational invariance of the correlation function, we calculate the spatially averaged correlation functions, i.e.,C t (t 1 , t 2 ) ≡ L −1 L 0 dx C ψ t (x; t 1 , t 2 ), which is equivalent to the corresponding correlation function above but helps to reduce the statistical error.
In Fig. S1, at linear system size L = 2 10 , the dependence of − log C t (t 1 , t 2 ) / C t (t 2 , t 2 ) on |t 1 − t 2 | for 8 different sets of parameters in the SCGLE are shown on a doublelogarithmic scale. The exponent β for different parameter choices is extracted from the slope of a selected portion of the corresponding curve of − log C t (t 1 , t 2 ) / C t (t 2 , t 2 ) on the double-logarithmic scale (cf. Fig. S1). For the noise levels σ = 20 −1 , 12 −1 , 10 −1 , 8 −1 , 5 −1 , the corresponding slopes are extracted by performing linear fits to the data points with |t 1 − t 2 | ∈ [10 3 , 10 4 ]. This gives rise to 2β = 0.61, for all the first four curves from below with σ = 20 −1 , 12 −1 , 10 −1 , 8 −1 , respectively, where we notice the typical KPZ scaling behavior, but no crossover effect to the anticipated exponential decay. This due to the fact that the crossover time scales corresponding to these noise levels are much larger than the time range available in our simulations.
For the fifth curve from below with σ = 5 −1 , one notices a considerably increased slope at the right end of the curve, clearly indicating the crossover effect which sets in within |t 1 − t 2 | ∈ [10 3 , 10 4 ]. This is reflected in the extracted exponent Figure S1.
(Color online) The dependence of − log C t (t 1 , t 2 ) / C t (t 2 , t 2 ) on |t 1 − t 2 | for 8 different sets of parameters in the SCGLE at linear system size L = 2 10  with 2β = 0.74. The anticipated exponential decay is expected to manifest itself at even larger |t 1 − t 2 |, which is however not accessible for our current computation resources.
For the first three curves from above with σ = 1, 2 −1 , 4 −1 , since the local gradient of the curve changing noticeably up to the last accessible data point with the largest value of |t 1 − t 2 |, denoted as ∆t max , the corresponding slopes are extracted by performing linear fits to the data points of a short segment of the curve with |t 1 − t 2 | ∈ [∆t max e −1/4 , ∆t max ] for σ = 4 −1 , 2 −1 , and an even shorter segment of the curve with |t 1 − t 2 | ∈ [∆t max e −1/5 , ∆t max ], for σ = 1. This gives rise to 2β = 1.0 for all the three curves, which clearly indicates the anticipated exponential decay. However, to further numerically verify this exponential decay in a larger time range is beyond our currently accessible computation resources.

PHASE RANDOM WALK MODEL
To see how the finite density of space-time vortices causes the exponential decay of the autocorrelation function, one can consider a simple "phase random walk" model in which the time evolution of the phase is assumed to be only determined by the phase jumps due to N uniformly distributed vortex cores at times nt v , where n = 1, . . . , N, with random charges W n = ±1 occurring with equal probability. The model can be formulated as where δΘ corresponds to the average amplitude of jumps of the phase field if a vortex core is crossed along the temporal direction. A straightforward calculation shows that ∆ θ (t 1 − t 2 ) = (δΘ) 2 |t 1 − t 2 | /t v , i.e., space-time vortices indeed lead to disordered behavior of the phase correlations beyond a time scale O(t v )

TECHNICAL DETAILS IN NUMERICAL SIMULATIONS
In numerical simulations, we discretize the SCGLE in space and time by ∆ x and ∆ t , respectively, and solve the discretized dynamical equation by using the semi-implicit algorithm developed in [3]. Each Gaussian white noise sequence ξ(x, t), which is used in the simulation for each corresponding stochastic trajectory of ψ(x, t), is generated by using the pseudorandom number generator employing the Mersenne Twister algorithm [4]. Typically, we choose ∆ x = 1, ∆ t ∈ [0.01, 0.1] in dimensionless units and the spatial size of the system L = 2 10 , which is large enough to avoid finite size effects in most simulations presented in the main text. However, for results corresponding to the low noise levels presented in the Fig. 2 in the main text, we have performed simulations on systems with spatial sizes up to 2 14 to avoid finite size effects that set in at long time. Typically, we simulate the SCGLE for a time period T ∈ [10 5 , 10 7 ] in dimensionless units, whose specific value depends on whether physical quantities of interest have reached their steady states. With the direct access to the time evolution of ψ(x, t) field for each stochastic trajectory simulated, the condensate temporal auto correlation function C ψ t (x; t 1 , t 2 ) ≡ ψ * (x, t 1 )ψ(x, t 2 ) and momentum distribution n q ≡ ψ * (q)ψ(q) can be calculated straightforwardly.
In order to calculate the temporal phase fluctuation function 2 , we need to further extract the dynamics of the phase field θ(x, t). At each time point t, we first assign θ(x, t) to be the complex argument field arg ψ(x, t) and then add a multiple of 2π to the value of θ(x, t) at each spatial point x in such a way that the phase difference between the new value of the phase at the current time point and the value of the phase at previous time step, i.e., θ(x, t) − θ(x, t − ∆ t ), assumes its value within the interval [−π, π). This procedure corresponds to a choice of the gauge for the phase field that eliminates the artificial discontinuous behavior associated with directly assigning θ(x, t) to be arg ψ(x, t). Afterwards, the temporal phase fluctuation function ∆ θ (t 1 − t 2 ) can be directly calculated.
To calculate the vortex density P v , for each stochastic trajectory, we first construct the discrete phase distribution θ i j ≡ (i∆x, j∆t), choosing ∆x = 1 and ∆t = 1, on the space-time plane of the size L × ∆T (with ∆T usually chosen to be 10 3 ). and then calculate the winding number of the θ i j field on each plaquette, where a vortex is identified whenever the winding number of the corresponding plaquette, i.e., the charge of the corresponding vortex, is non-zero. The vortex density P v is obtained by calculating N v /(L × ∆T ) , with N v being the total number of vortices on the space-time plane of the size L × ∆T for each stochastic trajectory. In Fig. S2, a typical snapshot of a vortex charge distribution on the space-time plane of the system in the vortex turbulence phase is shown, from which we notice the vortex charge distribution assumes an irregular structure.

THE FIRST ORDER TRANSITION LINE
In this section we present technical details concerning the investigation of the properties of the first order transition line. The dependence of P v on the nonequilibrium strengthλ at different noise levels with σ = 0.01, 0.011, 0.012, 0.013, 0.014, in the weak noise regime are shown in Fig. S3(a). At each noise level, we notice a sudden increase of P v by a few orders of magnitude whenλ is tuned across around 3.2, signifying a first order transition behavior. We estimate the critical valuẽ λ * for the transition at each noise level (cf. black dashed vertical line in Fig. S3(a)) by the arithmetic average of the values ofλ for the two data points right before and after the transition. This estimation gives rise toλ * = 3.23 with an error bar of ±0.07 for each noise level shown in Fig. S3(a), indicating that, within the numerical accuracy of our simulations,λ * is insensitive with respect to the noise level in the weak noise regime as shown in Fig. 3(c) in the main text.