Quentin Fontaine

Davide Squizzato

Florent Baboux

Ivan Amelio

Aristide Lemaître

Martina Morassi

Isabelle Sagnes

Luc Le Gratiet

Abdelmounaim Harouri

Michiel Wouters

I Carusotto

A Amo

M Richard

A Minguzzi

L Canet

S Ravets

J Bloch

dimensional polariton condensate

Revealing universal behaviors is a hallmark of statistical physics. Phenomena such as the stochastic growth of crystalline surfaces [START_REF] Krug | Universal finite-size effects in the rate of growth processes[END_REF], of interfaces in bacterial colonies [START_REF] Wakita | Self-affinity for the growing interface of bacterial colonies[END_REF], and spin transport in quantum magnets [START_REF] Ljubotina | Spin diffusion from an inhomogeneous quench in an integrable system[END_REF][START_REF] Ljubotina | Kardar-parisi-zhang physics in the quantum heisenberg magnet[END_REF][START_REF] Scheie | Detection of kardar-parisi-zhang hydrodynamics in a quantum heisenberg spin-1/2 chain[END_REF][START_REF] Wei | Quantum gas microscopy of kardar-parisi-zhang superdiffusion[END_REF] all belong to the same universality class, despite the great plurality of physical mechanisms they involve at the microscopic level. This universality stems from a common underlying effective dynamics governed by the non-linear stochastic Kardar-Parisi-Zhang (KPZ) equation [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF]. Recent theoretical works suggest that this dynamics also emerges in the phase of out-ofequilibrium systems displaying macroscopic spontaneous coherence [START_REF] Altman | Two-dimensional superfluidity of exciton polaritons requires strong anisotropy[END_REF][START_REF] Ji | Temporal coherence of one-dimensional nonequilibrium quantum fluids[END_REF][START_REF] He | Scaling properties of one-dimensional drivendissipative condensates[END_REF][START_REF] Zamora | Tuning across universalities with a driven open condensate[END_REF][START_REF] Comaron | Dynamical critical exponents in driven-dissipative quantum systems[END_REF][START_REF] Squizzato | Kardar-parisi-zhang universality in the phase distributions of one-dimensional exciton-polaritons[END_REF][START_REF] Amelio | Theory of the coherence of topological lasers[END_REF][START_REF] Ferrier | Searching for the kardar-parisi-zhang phase in microcavity polaritons[END_REF][START_REF] Deligiannis | Accessing kardar-parisi-zhang universality sub-classes with exciton polaritons[END_REF][START_REF] Mei | Spatiotemporal scaling of two-dimensional nonequilibrium exciton-polariton systems with weak interactions[END_REF]. Here, we experimentally demonstrate that the evolution of the phase in a driven-dissipative one-dimensional polariton condensate falls in the KPZ universality class. Our demonstration relies on a direct measurement of KPZ space-time scaling laws [START_REF] Family | Scaling of the active zone in the eden process on percolation networks and the ballistic deposition model[END_REF][START_REF] Halpin-Healy | Kinetic roughening phenomena, stochastic growth, directed polymers and all that. aspects of multidisciplinary statistical mechanics[END_REF], combined with a theoretical analysis that reveals other key signatures of this universality class. Our results highlight fundamental physical differences between out-ofequilibrium condensates and their equilibrium counterparts, and open a paradigm for exploring universal behaviors in driven open systems.

Universality is a powerful concept in statistical physics that allows describing critical phenomena on the basis of a few fundamental ingredients. At thermal equilibrium, models like the Ising model are pivotal in understanding the critical properties of a wide class of physical systems. Reversely, non-equilibrium systems lack a complete classification of their universal properties. In this context, the KPZ equation appears as a quintessential model to investigate non-equilibrium phenomena and phase transitions. Here, we provide the experimental demonstration that one-dimensional (1D) out-of-equilibrium condensates belong to the KPZ universality class.

The KPZ equation [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF] was originally proposed to describe the stochastic growth dynamics of an interface height h(r, t):

∂ t h(r, t) = ν∇ 2 h(r, t) + λ 2 [∇(h (r, t)] 2 + η(r, t). (1)
On the right, the first term corresponds to a smoothening diffusion, the second one to a nonlinear contribution leading to critical roughening, while the third one is a Gaussian white noise introducing stochasticity. The spatial and temporal correlation functions of h(r, t) exhibit power law behaviors, with critical exponents specific to the KPZ universality class [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF].

Currently available observations of KPZ dynamics have mainly focused on growing interfaces in classical systems [START_REF] Halpin-Healy | Kinetic roughening phenomena, stochastic growth, directed polymers and all that. aspects of multidisciplinary statistical mechanics[END_REF][START_REF] Krug | Origins of scale invariance in growth processes[END_REF][START_REF] Takeuchi | An appetizer to modern developments on the kardar-parisi-zhang universality class[END_REF]. Interestingly, KPZ physics is also under investigation in quantum magnets [START_REF] Ljubotina | Spin diffusion from an inhomogeneous quench in an integrable system[END_REF][START_REF] Ljubotina | Kardar-parisi-zhang physics in the quantum heisenberg magnet[END_REF][START_REF] Scheie | Detection of kardar-parisi-zhang hydrodynamics in a quantum heisenberg spin-1/2 chain[END_REF][START_REF] Wei | Quantum gas microscopy of kardar-parisi-zhang superdiffusion[END_REF].

Recent theoretical works predicted that the spatio-temporal evolution of the phase of a polariton condensate falls into the KPZ universality class [START_REF] Altman | Two-dimensional superfluidity of exciton polaritons requires strong anisotropy[END_REF][START_REF] Ji | Temporal coherence of one-dimensional nonequilibrium quantum fluids[END_REF][START_REF] He | Scaling properties of one-dimensional drivendissipative condensates[END_REF][START_REF] Zamora | Tuning across universalities with a driven open condensate[END_REF][START_REF] Comaron | Dynamical critical exponents in driven-dissipative quantum systems[END_REF][START_REF] Squizzato | Kardar-parisi-zhang universality in the phase distributions of one-dimensional exciton-polaritons[END_REF][START_REF] Amelio | Theory of the coherence of topological lasers[END_REF][START_REF] Ferrier | Searching for the kardar-parisi-zhang phase in microcavity polaritons[END_REF][START_REF] Deligiannis | Accessing kardar-parisi-zhang universality sub-classes with exciton polaritons[END_REF][START_REF] Mei | Spatiotemporal scaling of two-dimensional nonequilibrium exciton-polariton systems with weak interactions[END_REF]. However, unlike an actual interface height, the phase is defined periodically between 0 and 2π. This version of the KPZ equation is relevant for out-of-equilibrium systems developing macroscopic spontaneous coherence (lasers, arrays of coupled limit-cycles oscillators [START_REF] Lauter | From kardar-parisi-zhang scaling to explosive desynchronization in arrays of limit-cycle oscillators[END_REF]), and other systems like polar active smectic phases [START_REF] Chen | Universality for moving stripes: A hydrodynamic theory of polar active smectics[END_REF]. The compactness of the phase field results in a rich phase diagram comprising not only the KPZ phase but also other regimes characterized by the proliferation of topological defects [START_REF] Lauter | From kardar-parisi-zhang scaling to explosive desynchronization in arrays of limit-cycle oscillators[END_REF][START_REF] Chen | Universality for moving stripes: A hydrodynamic theory of polar active smectics[END_REF][START_REF] He | Space-time vortex driven crossover and vortex turbulence phase transition in one-dimensional driven open condensates[END_REF]. In this article, we explore experimentally the spatio-temporal dynamics of the first order coherence in a 1D polariton condensate. We observe the predicted coherence decay, and demonstrate the collapse of the data onto the universal KPZ scaling function. Our theoretical analysis evidences how the observed 1D KPZ physics is resilient to the presence of vortex-antivortex (V-AV) pairs.

Cavity polaritons are hybrid quasi-particles emerging in semiconductor cavities from the strong coupling between electronic excitations (excitons) in a quantum well and cavity photons [START_REF] Weisbuch | Observation of the coupled excitonphoton mode splitting in a semiconductor quantum microcavity[END_REF][START_REF] Carusotto | Quantum fluids of light[END_REF]. Owing to their excitonic fraction, polaritons can be created from a gain medium of incoherent excitons (the excitonic reservoir) through bosonic stimulated scattering. Due to photon leakage through the mirrors, the polariton dynamics is intrinsically out of equilibrium, its steady-state being the result of the balance between drive, relaxation and losses. Finally, polaritons can be laterally confined in lattices [START_REF] Schneider | Exciton-polariton trapping and potential landscape engineering[END_REF], enabling band structure engineering.

The 1D polariton condensates at play consist in the macroscopic occupation of a given state obtained by incoherently pumping the exciton reservoir [START_REF] Deng | Condensation of semiconductor microcavity exciton polaritons[END_REF][START_REF] Kasprzak | Bose-einstein condensation of exciton polaritons[END_REF]. Above a critical density, exciton stimulated scattering from the reservoir into this state triggers a spontaneous U (1) symmetry-breaking of the phase. The condensate and reservoir dynamics are described by two coupled equations [START_REF] Carusotto | Quantum fluids of light[END_REF]:

i ∂ ∂t ψ(x, t) = E( k) - i 2 γ( k) + g|ψ(x, t)| 2 + 2g R n R (x, t) + i 2 Rn R (x, t) ψ(x, t) + ξ(x, t) (2)
∂ ∂t n R (x, t) = P (x) -γ R + R|ψ(x, t)| 2 n R (x, t) (3)
Here, k = -i ∂/∂x is the momentum operator, ψ(x, t) = ρ(x, t) e iθ(x,t) the polariton condensate field with density ρ(x, t) and phase θ(x, t), E( k) the polariton dispersion, γ( k) the momentum-dependent decay rate, and g the polariton-polariton interaction strength. The exciton reservoir, with density n R (x, t), is pumped at rate P (x). Excitons either relax into the polariton condensate by stimulated scattering with rate R or decay following other channels at rate γ R . The term 2g R n R describes polariton repulsive interactions with reservoir excitons. It dominates the polariton blueshift close to threshold, and induces dephasing through inhomogeneous spectral broadening [START_REF] Love | Intrinsic decoherence mechanisms in the microcavity polariton condensate[END_REF]. Finally, ξ(x, t) describes Gaussian noise induced by drive and loss.

Ignoring interactions with the reservoir (g R = 0), previous theoretical studies have shown that the condensate phase θ(x, t) follows a KPZ equation [START_REF] Altman | Two-dimensional superfluidity of exciton polaritons requires strong anisotropy[END_REF][START_REF] Ji | Temporal coherence of one-dimensional nonequilibrium quantum fluids[END_REF][START_REF] He | Scaling properties of one-dimensional drivendissipative condensates[END_REF][START_REF] Grinstein | Temporally periodic phases and kinetic roughening[END_REF]. The condensate phase profile behaves as a classical interface (Fig. 1a), and develops KPZ spatio-temporal correlations characterized by the phase variance Var (∆θ(x, t)) = (∆θ(x, t)-∆θ(x, t) ) 2 (where . is the statistical averaging over different noise realizations, and ∆θ = θ(x, t)θ(x, t 0 ), t 0 being a reference time). Here, we derive the mapping to the KPZ equation for g R = 0 and obtain the KPZ parameters in terms of those entering Eqs. ( 2)-(3) (see [32]).

Experimentally probing KPZ correlations requires extended condensates to avoid finite size effects, a condition that was not fulfilled in early coherence measurements [START_REF] Roumpos | Power-law decay of the spatial correlation function in exciton-polariton condensates[END_REF][START_REF] Fischer | Spatial coherence properties of one dimensional exciton-polariton condensates[END_REF]. This requirement is demanding due to the development of a modulation instability, which fragments the condensate into mutually incoherent micron-sized puddles [START_REF] Bobrovska | Stability and spatial coherence of nonresonantly pumped exciton-polariton condensates[END_REF][START_REF] Daskalakis | Spatial coherence and stability in a disordered organic polariton condensate[END_REF][START_REF] Estrecho | Single-shot condensation of exciton polaritons and the hole burning effect[END_REF][START_REF] Bobrovska | Dynamical instability of a nonequilibrium exciton-polariton condensate[END_REF]. This instability originates from repulsive condensate-reservoir interactions, which result in effective attractive polariton-polariton interactions within the condensate and lead to its destabilization [START_REF] Smirnov | Dynamics and stability of dark solitons in exciton-polariton condensates[END_REF][START_REF] Liew | Instability-induced formation and nonequilibrium dynamics of phase defects in polariton condensates[END_REF].

A solution to tame this instability is to spatially separate the excitonic reservoir from the condensate [START_REF] Caputo | Topological order and thermal equilibrium in polariton condensates[END_REF], or to use negative mass polaritons, obtained by band engineering in a lattice [START_REF] Baboux | Unstable and stable regimes of polariton condensation[END_REF]. The negative mass changes the sign of the effective polariton-polariton interactions, thus restoring the condensate stability. We use this negative mass technique to generate stable 1D polariton condensates extending over distances as large as 100 µm (see Fig. 1b).

The sample consists in a semiconductor microcavity embedding quantum wells (Fig. 1c, and [32]). We use nanotechnology processes to fabricate 1D asymmetric Lieb lattices of coupled micropillars containing three sites per unit cell (see Fig. 1c, Methods Section, and [32]).

We incoherently populate the excitonic reservoir using a blue-detuned cw laser focused on a single lattice, with an elongated flat-top intensity profile.

The polariton emission analysed in momentum space below condensation threshold (see Fig. 2a) evidences the lattice band structure emerging from the hybridization of the discrete modes confined in each micropillar. Above a threshold P th = 50 mW, the emission becomes peaked at the top of a negative mass band (Fig. 2b). This feature, together with the non-linear increase of the emission intensity (Fig. 2c), indicates the onset of polariton condensation. The condensate emission intensity in real space at P = 1.1 P th reveals an extended and regular intensity profile envelope, as expected in absence of modulation instability (Fig. 1b). Although we use a flat-top intensity profile for the excitation beam, the condensate profile is rather peaked, which hints toward localization induced by disorder in the lattice.

We define the first order correlation evaluated between points separated in space by ∆x = 2x and delayed by ∆t:

g (1) (∆x, ∆t) = ψ * (x, t 0 )ψ(-x, t 0 + ∆t) |ψ(x, t 0 )| 2 |ψ(-x, t 0 + ∆t)| 2 . ( 4
)
Neglecting density-density and density-phase correlations, we show that g (1) (∆x, ∆t) ≈ exp {-Var [∆θ(∆x, ∆t)] /2} (see [32]). We thus expect KPZ universal scaling to show up as stretched exponentials in the coherence decay: |g (1) (∆x, ∆t 0 )| ∝ exp [-(∆x/λ) 2χ /2] (for fixed ∆t 0 ) and |g (1) (∆x 0 , ∆t)| ∝ exp -(∆t/τ ) 2β /2 (for fixed ∆x 0 ), where χ and β are the universal KPZ critical exponents while λ and τ are two non-universal parameters. In 1D, the "roughness" exponent χ is equal to 1/2 and the "growth" exponent β to 1/3 [START_REF] Family | Scaling of the active zone in the eden process on percolation networks and the ballistic deposition model[END_REF][START_REF] Halpin-Healy | Kinetic roughening phenomena, stochastic growth, directed polymers and all that. aspects of multidisciplinary statistical mechanics[END_REF].

While χ = 1/2 is common to several universality classes for 1D systems (such as linear systems described with Bogoliubov theory [START_REF] Edwards | The surface statistics of a granular aggregate[END_REF][START_REF] Wouters | Absence of long-range coherence in the parametric emission of photonic wires[END_REF]), β = 1/3 is an unambiguous signature of KPZ physics.

The condensate coherence |g (1) | is measured using Michelson interferometry (Fig. 2d). The field emitted by the condensate at time t 0 , E(x, t 0 ) ∝ Ψ(x, t 0 ), and the one emitted at the mirror-symmetric point at time t 0 + ∆t, E(-x, t 0 + ∆t) ∝ Ψ(-x, t 0 + ∆t), are overlapped (∆t is the delay in the two interferometer arms). The resulting intensity pattern exhibits well-contrasted interference fringes over the whole condensate (Fig. 2e). This indicates the emergence of extended spatial coherence as the fringe contrast gives a direct visualization of the degree of coherence between fields emitted at two points spatially separated by ∆x = 2x and delayed by ∆t. More specifically, |g (1) (∆x, ∆t)| is determined from the fringe visibility V (∆x, ∆t) and from the intensity distributions |E(x, t 0 )| 2 and |E(-x, t 0 + ∆t)| 2 measured separately, using

2|g (1) (∆x, ∆t)| |E(x, t 0 )| 2 |E(-x, t 0 + ∆t)| 2 = V (∆x, ∆t) |E(x, t 0 )| 2 + |E(-x, t 0 + ∆t)| 2 .
(

) 5
The result is shown in Fig. 2f.

We first focus on the temporal decay of the coherence. To search for the growth exponent, we compute the temporal derivative D t = -2 ∂ log |g (1) (0, ∆t)| /∂∆t from our dataset. According to KPZ theory, this derivative scales as a power law with exponent 2β -1 = -1/3. In the inset of Fig. 3a, we identify such scaling throughout the temporal window 15 ps ≤ ∆t ≤ 80 ps. Equivalently, we observe in the main panel a linear increase of -2 log |g (1) (0, ∆t)| as a function of ∆t 2/3 over the same window (grey shaded area), which demonstrates a key feature of KPZ dynamics. At short timescales, the deviation from KPZ scaling and the saturation of |g (1) | are due to an incoherent background (spectrally broad photoluminescence from uncondensed states) that hides the onset of KPZ fluctuations. For ∆t ≥ 80 ps, -2 log |g (1) (0, ∆t)| also departs from KPZ scaling and follows a super-linear behavior that we attribute to slow reservoir population fluctuations [START_REF] Love | Intrinsic decoherence mechanisms in the microcavity polariton condensate[END_REF]. We now perform a similar analysis in the spatial domain. The results are shown in Fig. 3b. The spatial derivative D x = -2 ∂ log |g (1) (∆x, 0)| /∂∆x exhibits a plateau within the spatial window 30 µm ≤ ∆x ≤ 60 µm, in agreement with χ = 1/2 (inset). The roughness exponent χ = 1/2 also shows up in the linear increase of -2 log |g (1) (∆x, 0)| as a function of ∆x, over the same window (grey shaded area in the main panel). When approaching condensate edges (∆x ≥ 60 µm), the coherence decays faster as a consequence of enhanced fluctuations at smaller polariton density (see [START_REF] Chiocchetta | Non-equilibrium quasi-condensates in reduced dimensions[END_REF]). Pushing further this data analysis, we fit the coherence decay curves with stretched exponentials and deduce experimental values for the scaling exponents: χ exp = 0.51 ± 0.08 and β exp = 0.36 ± 0.11. The uncertainty on β exp allows us to discriminate between the different universality classes relevant for our system, as the KPZ value β = 1/3 remains the only one lying within the 95% confidence interval on β exp (see [32]).

We now search for KPZ signatures over the whole space-time correlation map. We select all data points where 0.57 ≤ |g (1) | ≤ 0.75, the range where we evidence KPZ scaling at ∆t = 0. This space-time window is shown in the bottom inset of Fig. 3c. For this subset of data points, we plot in Fig. 3c the value of -2 log κ|g (1) (∆x, ∆t)| /∆t 2β as a function of the rescaled coordinate y = ∆x/∆t 1/Z , where Z = β/χ = 2/3 and κ is a normalization factor (see [32]). Strikingly, all these data points collapse onto the scaling function

F = C 0 F KPZ (y/y 0 )
, where F KPZ is the tabulated dimensionless KPZ universal scaling function [START_REF] Prähofer | Exact scaling functions for one-dimensional stationary kpz growth[END_REF][START_REF] Spohn | Nonlinear fluctuating hydrodynamics for anharmonic chains[END_REF]. We use the non-universal constants C 0 and y 0 as fitting parameters in order to vertically and horizontally shift F onto the collapsed data points. This result demonstrates that 1D polariton condensates indeed belong to the KPZ universality class. In order to reinforce the generality of this conclusion, we performed the same measurement and analysis on a different 1D lattice with four sites per unit cell. We also found a KPZ space-time window where all data points collapse onto the universal scaling curve [32]. To complete the picture, we carried out the same analysis at higher excitation powers, and found for both lattices that the spatio-temporal KPZ window shrinks for increasing P/P th and eventually disappears when P/P th > 1.2 [32].

With the prospect of confronting these experimental data to simulations, we compute the phase evolution of the polariton condensate by numerically solving Eqs. ( 2)- [START_REF] Ljubotina | Spin diffusion from an inhomogeneous quench in an integrable system[END_REF]. We obtain some microscopic parameters from the experiment, namely: the dispersion relation E(k), the polariton blueshift 2g R n R , the linewidth γ(k), and the threshold-normalized laser power P/P th . For the reservoir relaxation rate R and decay rate γ R , we choose values within some realistic range yielding the best agreement with the measured |g (1) |. We also adjust the pump size to make the condensate intensity profile |Ψ(x)| 2 comparable with the experimental one. Finally, we neglect the polariton-polariton interaction energy g|Ψ| 2 and take g = 0 in all simulations. As such, this model also applies to spatially extended lasers in the weak coupling regime. The calculated |g (1) num | data are reported in Fig. 3a-b, showing excellent agreement with the experiment. The short and long time behaviour is reproduced, together with the shrinking of the KPZ window when the excitation power is increased (see [32]).

We then perform on the numerical data the same analysis as on the experimental ones.

We plot -2 log κ|g (1) num (∆x, ∆t)| /∆t 2β as a function of y, selecting the points for which 0.57 ≤ |g (1) num | ≤ 0.75. The result is shown in Fig. 3c (top inset), together with the KPZ scaling function F = C 0 F KPZ (y/y 0 ), using for C 0 and y 0 the same values as for the experimental data. The simulated data align to the scaling function, thus fully validating our model.

To deepen our insight into the phase dynamics, we now analyze its stochastic behavior in the numerical simulations. Fig. 4a shows an example of a phase map θ(x, ∆t) corresponding to a given realization of the noise (others are shown in [32]). For better visualization, we plot the wrapped phase after subtracting the dynamical phase ω 0 ∆t ( ω 0 being the condensate energy). We observe two kinds of phase variations: small amplitude fluctuations and fast (scarce) phase jumps. These jumps are associated to pairs of close-by spatio-temporal vortices with opposite circulation (see inset), that we name V-AV pairs. To analyse the effect of these V-AV pairs on the phase dynamics, we show in Fig. 4b the unwrapped phase temporal evolution at x = 0 (horizontal line in Fig. 4a). The phase evolution exhibits plateaus with small amplitude phase fluctuations, separated by phase jumps of approximately 2π, occurring on a fast timescale (∼ 1 ps) when passing through a V-AV pair. Note that for the regime of parameters explored here, almost all vortices appear in V-AV pairs. For higher noise or stronger interactions, activation of single vortices is expected and would lead to other dynamical regimes [START_REF] He | Space-time vortex driven crossover and vortex turbulence phase transition in one-dimensional driven open condensates[END_REF].

We now focus on showing that small amplitude phase fluctuations follow KPZ scaling laws by computing the phase variance Var [∆θ], and that the presence of V-AV only weakly affects |g (1) num |. Since Var [∆θ] is extremely sensitive to phase jumps, we select for each trajectory a 100 ps wide window where the phase undergoes the smallest amount of jumps. When few phase jumps remain in the selected window, we filter them out by adding at every vortex location an other one of opposite charge (grey line in Fig. 4b). Note that we discard 5% of the realizations, where vortices proliferate (see [32]). We then compute the phase variance Var [∆θ VF ] over the set of vortex-free (VF) time windows. The result is plotted in Fig. 4c together with the values of -2log |g (1) num | (computed without filtering the V-AV pairs). Both quantities exhibit the KPZ power-law scaling over the same time window (grey-shaded area), as further illustrated on Fig. 4c (inset) where we plot their time derivative. This result reveals that the first order coherence is indeed a good observable to probe the KPZ dynamics of the condensate phase, even in presence of V-AV pairs.

Another striking signature of KPZ physics lies in the fact that phase fluctuations are governed by a probability distribution, whichunlike the normal distributionis skewed and exhibits markedly different tails. In Fig. 4d, we show the calculated probability distribution of the unwrapped ∆θ(0, ∆t), computed over all trajectories (thus including vortices) for ∆t = 50 ps, i.e. in the center of the KPZ window. All trajectories that have not crossed any V-AV pair contribute to the first peak in the distribution. The second peak corresponds to trajectories which have crossed one V-AV pair before reaching ∆t = 50 ps and have thus undergone one phase jump close to 2π. Strikingly the first two peaks are skewed and well reproduced by the Tracy-Widom Gaussian Orthogonal Ensemble (GOE) distribution associated to the flat subclass (see [32]). Cumulating data for various ∆t, we obtain an agreement with the Tracy-Widom GOE distribution over six decades (see [32] for details on the analysis). The third peak corresponds to few realizations displaying two phase jumps.

The lack of statistics prevents precise analysis of its shape. Our simulations highlight that V-AV pairs only modify the probability distribution by adding replicas of the main peak without significantly changing their shape. Moreover, they confirm that in the regime where a low density of V-AV pairs stochastically shows up, KPZ dynamics is not destroyed but occurs piece-wise in the spatio-temporal domain.

To conclude, both our experimental and theoretical analysis prove that KPZ scaling laws are present in the decay of the first-order coherence of 1D driven-dissipative polariton condensates. Our findings apply to any spatially extended driven open systems subject to gain and loss and characterized by a U (1) symmetry breaking. Our work opens many new challenges to be addressed in the future. In 1D, while our results highlight the striking resilience of KPZ physics to V-AV pairs, regimes at higher noise strength or higher non-linearity remains to be explored [START_REF] He | Space-time vortex driven crossover and vortex turbulence phase transition in one-dimensional driven open condensates[END_REF]. Investigating different KPZ universality sub-classes, predicted for various geometries of the condensate environment [START_REF] Deligiannis | Accessing kardar-parisi-zhang universality sub-classes with exciton polaritons[END_REF], is now also within reach, when engineering the geometry of the condensate environment. Beyond 1D, exciton-polariton lattices offer exciting perspectives for the exploration of KPZ physics in 2D, where an experimental realization is highly sought-after [START_REF] Zamora | Tuning across universalities with a driven open condensate[END_REF][START_REF] Comaron | Dynamical critical exponents in driven-dissipative quantum systems[END_REF][START_REF] Squizzato | Kardar-parisi-zhang universality in the phase distributions of one-dimensional exciton-polaritons[END_REF][START_REF] Amelio | Theory of the coherence of topological lasers[END_REF][START_REF] Ferrier | Searching for the kardar-parisi-zhang phase in microcavity polaritons[END_REF][START_REF] Deligiannis | Accessing kardar-parisi-zhang universality sub-classes with exciton polaritons[END_REF], and the role of topological defects still actively debated [START_REF] Zamora | Tuning across universalities with a driven open condensate[END_REF]. An experimental implementation involving polariton condensates would enable to test the different models and serve as a general analog simulator of complex physical systems belonging to KPZ universality class.

Method

The We incoherently populate the excitonic reservoir using a non-resonant continuous-wave laser at 740 nm (reflectivity minimum of the Bragg mirror). A spatial light modulator (SLM) enables shaping the excitation spot into a 125 µm long flat-top beam in the lattice direction.

Note that our experiment is performed under truly continuous-wave excitation conditions, that is, without any chopper. The polariton emission leaking out through the cavity top mirror is analyzed in space, momentum (along the lattice direction x) and frequency with a monochromator coupled to a CCD camera.

We retrieve the condensate first order coherence using Michelson intererometry. A 2-mirror retro-reflector mounted on a step-motorized translation stage in one of the interferometer arms enables us to overlap on a CCD camera the field E(x, t 0 ) emitted by the condensate at time t 0 and position x, with E(-x, t 0 +∆t), the field emitted at t 0 +∆t and position -x (∆t is the delay introduced between the interferometer arms by translating the retro-reflector). In order to probe the temporal scaling of the condensate coherence, we scan the retro-reflector position over ∆L = 5 cm, corresponding to a maximum time delay of ∆t = 2∆L/c = 330 ps.

During such a scan, we set the camera exposure time to 1s and acquire a series of 250 images. (1) |, retrieved from the fringe visibility using Eq. ( 5), as a function of ∆x and ∆t. Top inset: Same as in the main panel but for numerical data. On both graphs, the black line corresponds to the KPZ scaling function F = C 0 F KPZ (y/y 0 ), adjusted to the experimental data by tuning the values of C 0 and y 0 . Experimental and simulated |g (1) | datasets are normalized using κ = 1.13. In all panels, the numerical data are averaged over 10 4 realizations of the noise.

I. OVERVIEW

In this Supplemental Material, we provide additional information on the experiments and on the numerical simulations, as well as additional discussion of the results. In Sec. II, we analytically derive the mapping from our two-coupled equation model for the dynamics of the condensate field and reservoir density to the KPZ equation for the phase dynamics.

We precisely relate the g (1) first-order correlation function to the phase-phase correlations.

In Sec. III, we provide all the information on the experimental set-up and measurements, and we report complementary experimental results obtained on a symmetric Lieb lattice. In Sec. V we perform an in-depth analysis of the phase dynamics and of the effect of space-time vortices.

Beyond all the necessary discussion, let us emphasize below the main results reported in this material:

• We establish the mapping to the KPZ equation for a more general and realistic model than previous studies in (Sec. II A).

• We consolidate the validity of our experimental findings by reproducing them in a different lattice featuring condensation in a different type of bands (Sec. III E).

• We demonstrate that the measured scaling behavior of g (1) directly reflects the KPZ scaling of the phase (Sec. V B),

• We analyze the effects of space-time vortices, and explain why KPZ dynamics can be resilient to their presence (Sec. V C).

II. THE THEORETICAL MODEL: EMERGENCE OF KPZ DYNAMICS IN IN-COHERENTLY PUMPED POLARITONS

A. The driven-dissipative Gross-Pitaevskii equation under incoherent pumping

We consider an out-of-equilibrium polariton condensate created in a one-dimensional lattice. Since the relevant dynamics occurs at low energy, we restrict the description to an effective single-band model, neglecting the contribution of the other lattice bands. We describe the polariton condensate wavefunction by the classical field ψ(x, t) at position x and time t. The excitation of the polariton condensate is modeled by introducing an external pump P (x) filling an incoherent excitonic reservoir of density n R (x, t). The reservoir excitons either relax into the polariton condensate by stimulated scattering with rate R or decay via other channels with total rate γ R [START_REF] Porras | Polariton dynamics and bose-einstein condensation in semiconductor microcavities[END_REF]. We describe the polariton-polariton and excitonpolariton interactions as contact interactions of strength g and g R respectively. The coupled equations for the condensate and reservoir read [START_REF] Wouters | Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons[END_REF]:

     i ∂ t ψ = E( k) - i 2 γ( k) + g|ψ| 2 + 2g R n R + i 2 R n R ψ + ξ ∂ t n R = P -(γ R + R|ψ| 2 )n R , (6)
where the first equation is a generalized, stochastic Gross-Pitaevskii equation (gGPE) for the condensate and the second one is a rate equation for the reservoir. In Eq. ( 6) k = -i ∂/∂x is the momentum operator and ξ is a white noise with correlations ξ(x, t)ξ * (x , t ) = 2ξ 0 δ(xx )δ(tt ). In the vicinity of k = 0, we approximate the lattice dispersion by

a parabola E(k) = E 0 + 2 k 2 /2m
. The polariton linewidth γ(k) is taken as momentum dependent, consistently with the experimental observations (see Sec. III C below) and is well approximated by γ(k 0) γ 0 + γ 2 k 2 near k = 0. The momentum-dependent linewidth plays a very important role in our model as it ensures the stability of the polariton condensate in our simulations [START_REF] Baboux | Unstable and stable regimes of polariton condensation[END_REF]. Note that it also has a crucial impact on the edge dynamics of topological lasers, where k-dependent losses naturally occur from the k-dependent confinement of the edge mode [START_REF] Amelio | Theory of the coherence of topological lasers[END_REF][START_REF] Loirette-Pelous | Linearized theory of the fluctuation dynamics in two-dimensional topological lasers[END_REF]. For the variance of the noise used in Eq. ( 6), we take the value ξ 0 = R 2 n R , that represents quantum noise due to pumping within the truncated Wigner picture. In this way, the correlators of the quantum field can be extracted from the noise-averaged spatio-temporal correlators of ψ(x, t) [START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Wouters | Stochastic classical field model for polariton condensates[END_REF].

B. Mapping to the KPZ equation

In the following, we use the density-phase representation of ψ(x, t) within the rotating frame of the condensate: ψ(x, t) = ρ(x, t) exp i(θ(x, t)ω 0 t) , where

ω 0 = g R R γ 0 2 + g g R γ R
γ 0 (p-1) is the condensate emission frequency. We focus on the dynamics of the system at small momentum and frequency. For the sake of simplicity, we therefore perform our analysis using the parabolic approximations of both E(k) and γ(k), introduced in section II A. In terms of the phase and density fields, the Laplacian and time-derivative operators read:

∂ t ψ =ψ 1 2 ρ -1 ∂ t ρ + i∂ t θ -iω 0 , (7)
∇ 2 ψ(x, t) = ψ(x, t) - 1 4 ρ(x, t) -2 (∇ρ(x, t)) 2 + 1 2 ρ(x, t) -1 ∇ 2 ρ(x, t) -(∇θ(x, t)) 2 + iρ(x, t) -1 ∇ρ(x, t) • ∇θ(x, t) + i∇ 2 θ(x, t) ≡ ψ(x, t) D[ρ, θ] , (8)
where in the last line we have introduced the differential operator D[ρ, θ] for simplicity in the notation. Equations (6) turn into a set of three coupled equations for the real-valued fields θ(x, t), ρ(x, t) and n R (x, t),

             ∂ t θ = ω 0 + 2m {D[ρ, θ]} + γ 2 2 {D[ρ, θ]} - g ρ -2 g R n R + { ξ} 1 2ρ ∂ t ρ = - 2m {D[ρ, θ]} + γ 2 2 {D[ρ, θ]} + 1 2 R n R -γ 0 + { ξ} ∂ t n R = P -(γ R + R ρ) n R (9)
where {•}, {•} stand for the real-and imaginary-part and ξ = -ie -iθ+iω 0 t ρ -1/2 ξ.

The existence of a Goldstone mode entails that phase fluctuations dominate in the long time, large-distance regime. Conversely, the polariton and reservoir densities are subject to a restoring force that makes them relax, within a short timescale, to their stationary values ρ 0 , n R,0 . It is thus convenient to perform the following decomposition: δρ = ρ(x, t)ρ 0 and δn R = n R (x, t)n R,0 . We now assume that δρ/ρ 0 1 and δn R /n R,0 1, and that these fluctuations are stationary, i.e. we neglect ∂ t δρ and ∂ t δn R . This is justified when the time scales of the reservoir and condensate density fluctuations are well separated from the ones of the phase fluctuations. This is similar in spirit to the usual decoupling approximation [START_REF] Altman | Two-dimensional superfluidity of exciton polaritons requires strong anisotropy[END_REF],

but on the two equations for the densities. We also neglect the spatial dependence of the condensate density fluctuations, which leads to the following set of equations, describing the effective dynamics of the system:

               ∂ t θ = - 2m (∇θ) 2 + γ 2 2 ∇ 2 θ - g δρ -2 g R δn R + { ξ} δn R = 2 R 2m ∇ 2 θ + γ 2 2 (∇θ) 2 -{ ξ} δρ = -p γ R γ 0 δn R . (10)
After a simple substitution, the equation governing the phase evolution becomes:

∂ t θ = γ 2 2 -u g R R m ∇ 2 θ - 2m + u g R R γ 2 (∇θ) 2 + η ≡ ν∇ 2 θ + λ 2 (∇θ) 2 + η (11)
with

η(x, t)η(x , t ) = ξ 0 ρ 0 1 + 4 u g R R 2 δ(x -x )δ(t -t ) ≡ 2Dδ(x -x )δ(t -t ) , (12)
and

u = 2 -p g g R γ R γ 0 . (13)
Equation ( 11) is the KPZ equation for the phase. The terms neglected during the derivation of this equation may slightly renormalize the KPZ parameters, but they do not drive the phase dynamics out of the KPZ universality class. The numerical results (see Sec. V below)

show that their effect is negligible in our case.

Equations ( 11) and ( 12) allow us to obtain the expression of the KPZ parameters ν, λ and D in terms of the microscopic parameters entering Eq. ( 6):

ν = γ 2 2 -u g R R m , λ = -2 2m + u g R R γ 2 , D = R n R,0 4ρ 0 1 + 4 u g R R 2 , (14)
with u given by Eq. ( 13). Note that the effective diffusivity ν must be positive in order for the linear limit of the KPZ equation, i.e. the Edwards-Wilkinson equation (EW), to be stable. Therefore, the expression of ν provides important insights into the role played by the non-trivial k-dependence of the linewidth, expressed by the parameter γ 2 , and by the sign of the mass. For the parameters used in our simulations, u is positive (see section IV).

The negative sign of the polariton mass is thus crucial in stabilizing the system [START_REF] Baboux | Unstable and stable regimes of polariton condensation[END_REF]. In the case of a positive mass, a non-zero value of γ 2 is absolutely necessary in order for the phase to be stable: a vanishing value of γ 2 would yield a negative value for ν, hence an instability of the KPZ equation.

Finally, we would like to stress that our derivation of the mapping between the generalized Gross-Pitaevskii equation for the polariton condensate and the KPZ equation for the phase is more general and realistic than those found in previous studies [8-13, 15-17, 52]: we formulate the decoupling approximation for the three-equation system ( 9) and we include in the gGPE the reservoir-induced blue-shift term 2g R n R . Both aspects are crucial to provide a faithful description of the experiment.

C. Connection between the condensate first-order correlation and the two-point phase-phase correlations

In this section, we detail the link between the condensate first-order correlation function g (1) , which is measured experimentally, and the phase-phase correlation function, which displays universal spatio-temporal KPZ scaling. In particular, we derive the conditions required in order to ensure that the scaling behavior of g (1) reflects the underlying KPZ dynamics of the condensate phase.

The general definition of the first-order correlation reads:

g (1) (∆x, ∆t) = ψ * (x, t 0 )ψ(-x, t 0 + ∆t) ρ(x, t 0 ) ρ(-x, t 0 + ∆t) = ρ(x, t 0 )ρ(-x, t 0 + ∆t)e i∆θ(∆x,∆t) ρ(x, t 0 ) ρ(-x, t 0 + ∆t) . ( 15
)
In the left-hand side, we omitted the dependence on t 0 due to the stationarity of the condensate dynamics. In our work, we are interested in accessing the phase dynamics from the field-field correlator. If we assume that the dynamics of the phase is decoupled from the one of the density, we get:

g (1) (∆x, ∆t) = ρ(x, t 0 )ρ(-x, t 0 + ∆t) ρ(x, t 0 ) ρ(-x, t 0 + ∆t) e i∆θ(∆x,∆t) . ( 16
)
We then decompose the density field into a mean-field and a fluctuating contribution, ρ(x, t) = ρ 0 + δρ(x, t), and assume that δρ(x, t)/ρ 0 1. We expand both the numerator and the denominator in the right-hand side of Eq. ( 16), which become, to linear order in δρ/ρ 0 :

ρ(x, t 0 )ρ(-x, t 0 + ∆t) ρ 0 1 + 1 2 δρ(x, t 0 ) ρ 0 1 + 1 2 δρ(-x, t 0 + ∆t) ρ 0 ρ 0 + 1 2 ( δρ(x, t 0 ) + δρ(-x, t 0 + ∆t) ) , (17)
and ρ(x, t 0 ) ρ(-x, t 0 + ∆t)

ρ 0 1 + 1 2 δρ(x, t 0 ) ρ 0 1 + 1 2 δρ(-x, t 0 + ∆t) ρ 0 ρ 0 + 1 2 ( δρ(x, t 0 ) + δρ(-x, t 0 + ∆t) ) . (18)
In this limit, the density terms in Eq. ( 16) simplify and one hence gets

g (1) (∆x, ∆t) = exp [i∆θ(∆x, ∆t)] . (19)
Furthermore, for small fluctuations of the phase, we can use the cumulant expansion to get:

|g (1) (∆x, ∆t)| 2 exp(-∆θ(∆x, ∆t) 2 + ∆θ(∆x, ∆t) 2 ) ≡ exp(-Var [∆θ(∆x, ∆t)]) , (20)
and hence

-2 log(|g (1) (∆x, ∆t)|) Var [∆θ(∆x, ∆t)] . (21)
It is instructive to study the validity of approximation [START_REF] Krug | Origins of scale invariance in growth processes[END_REF]. At O(∆θ 4 ) we have |g (1) (∆x, ∆t)| 2 = 1 -∆θ(∆x, ∆t)

We thus expect the two quantities to differ significantly when O(∆θ 4 ) fluctuations become comparable to 1 -Var [∆θ]. The effect of the density-density and density-phase correlations are studied in Sec. V B. We show that they do not affect the KPZ scaling for all time delays within the time window where KPZ scaling is observed, thus supporting the corresponding assumption in the derivation of the relation [START_REF] Halpin-Healy | Kinetic roughening phenomena, stochastic growth, directed polymers and all that. aspects of multidisciplinary statistical mechanics[END_REF]. The validity of Eq. ( 21) is also discussed in the same section.

B. From a single micropillar to lattices

Micropillars constitute the elementary building block of the lattices we use. In such a structure, the electric field is confined in all directions: longitudinally by the cavity mirrors, transversely by the large refractive index mismatch between AlGaAs and vacuum. In the transverse plane, polaritons are thus confined through their photonic component in a quasiinfinite circular potential. This confinement yields discrete energy modes whose spatial shape are similar to the hydrogen atomic orbitals. The lowest energy mode has a single bright lobe and thus corresponds to a S-state; the next two modes correspond to P-states; and so on.

As mentioned in section III A, the unit cell of a 1D asymmetric Lieb lattice contains three sites (labeled A, B and C in Fig. 5a), linked by the coupling constant t. In the quasicontinuum limit where several unit cells are arranged along a 1D lattice, this coupling between sites yields the hybridization of the pillar S-orbitals into three dispersive S-bands, gapped one from the other. The same reasoning enables describing the appearance of six higher energy P-bands, resulting from the hybridization of the pillar P-orbitals.

The 1D symmetric Lieb lattice contains four sites (labeled A, B, C and D in Fig. 5b), and thus presents four dispersive S-bands, and eight dispersive P-bands.

C. Asymmetric lattice characterization -Microscopic parameters

Low-power photoluminescence spectrum

Linear spectroscopy enables visualizing the band structure of our lattices, from which we can extract some of the parameters entering our numerical simulations. The inset of Fig. 6a shows the far-field emission (in TM polarization, parallel to the lattice axis) of the asymmetric Lieb lattice at low excitation power (P/P th ≈ 0.5). The bottom three bands visible on this image correspond to the three lattice S-bands. Above, we also see the first P-band, separated by a small gap from the upper S-band, at the top of which condensation takes place. Fitting the latter with a Lorentzian lineshape for all wave-vectors lying in the first Brillouin zone enables us to retrieve all at once:

-the polariton dispersion E(k) (see Fig. 6a), from which we extract the polariton mass m = -3.3 × 10 -6 m e (where m e is the electron mass);

-the polariton group velocity v g (k) (see Fig. 6b), obtained by differentiating the dispersion;

-the spectral linewidth γ spe (k) (light blue points in Fig. 6c), from which we obtain an estimate of the polariton linewidth at k = 0: γ spe (0) ≈ 80 µeV.

The k = 0 value of the measured spectral linewidth appears to be relatively large compared to the 22 µeV nominal linewidth expected for this structure. This is most probably due to electrostatic fluctuations in the sample during the integration time (∼ 60 s), which induce a spectral wandering of the emission energy through the polariton excitonic component [START_REF] Kuhlmann | Charge noise and spin noise in a semiconductor quantum device[END_REF].

This leads, in turn, to an inhomogeneous broadening of the polariton linewidth.

Propagation measurement

We can get a better estimate of the polariton linewidth γ(k) by probing in real space the energy resolved propagation of polaritons along the lattice, under localized excitation.

Depending on their wave-vector, polaritons travel away from the excitation spot, with a constant group velocity |v g (k)|. Because of their finite lifetime, this propagation results in an exponential decrease of the photoluminescence intensity along the lattice direction (I(x) ∝ exp{-|x|/L x }), as shown in Fig. 7a. Fitting this decay at different energies E(k)

allows us to retrieve the polariton linewidth γ pro = v g (k)/L x as function of k (red dots in Fig. 7b).

This method has the advantage of being less sensitive to charge fluctuations. Indeed, polaritons leave the pumping area with a given initial energy E 0 , setting the group velocity at which they travel. This group velocity remains constant over the whole propagation as charge fluctuations (i) mainly affect the reservoir energy locally (under the pump spot) and

(ii) occur on a time scale much larger than the polariton lifetime. Therefore, the propagation length L x only depends on E 0 regardless of the exciton energy at the time at which the polariton was emitted. Consequently, we assume that the measurement of γ pro is less affected by the wandering of the exciton energy and, moreover, that it almost corresponds to the Lorentzian contribution to γ spe .

The Gaussian contribution γ g to the spectral linewidtharising from the inhomogeneous broadeningcan then be retrieved using the following approximation [START_REF] Olivero | Empirical fits to the voigt line width: A brief review[END_REF]: We notice that γ g (diamonds in Fig. 7b) is nearly constant and equal to 73 µeV in the range 0.1 π/a < k < 0.7 π/a where the propagation measurement is reliable (outside of this range, v g is too small to properly extract γ pro ). Assuming that γ g remains constant over the Brillouin zone, we can finally remove the contribution of the inhomogeneous broadening to the spectral linewidth and obtain a better estimate of the polariton linewidth γ (red stars in Fig. 7c). We find γ(0) = 40 ± 10 µeV, which is reasonable compared to the nominal linewidth given earlier. The large errorbars on the data points (red shaded area) mainly come from the uncertainty on the measurement of γ g . We also show in Fig. 7d a comparison between the experimental data and the fit of the linewidth behavior used in the numerical simulations (see section IV), and its parabolic approximation around k = 0 (γ = γ 0 + γ 2 k 2 ), used in the derivation of the mapping in Sec. II B (red dotted line). The optimal parameters found in our simulations are given by γ 0,th = 48.5 µeV and γ 2,th = 1.6 × 10 4 µeV.µm 2 .

γ spe = 0.535 γ pro + 0.217 γ 2 pro + γ 2 g . ( 24
)
D. Optical setup and data analysis This arrangement of the interferometer enables us to tune the interfringe spacing of the resulting interference pattern, as it allows to control the incident wave-vector of both beams before the last lens as well as their relative spacing. In order to probe the temporal scaling of the condensate first-order correlation function, we typically scan the retroreflector position over a distance of ∆L = 5 cm, corresponding to a maximum time delay of ∆t = 2∆L/c = 330 ps. During such a scan, we set the camera exposure time to 1 s and acquire a serie of 250 images. The zero delay position (∆L = 0) has been calibrated beforehand by sending white light through the interferometer.

Data analysis procedure

In our experimental setup, the condensate image (reference arm) is overlapped with image at the mirror-symmetric point with respect to a plane orthogonal to the lattice (retroreflector arm). The resulting interference pattern (at ∆t = 0) is shown in Fig. 8a and Fig. 12d for the asymmetric and symmetric Lieb lattices respectively. At each point r = (x, y) of the image plane, the intensity I c (x, y) is given by the interference between the fields emitted at

x and -x in the sample plane. Dropping the y coordinate, we thus expect that:

I c (x, ∆t) = 1 4 I(x) + I(-x) + 2 I(x)I(-x) |g (1) (∆x, ∆t)| cos (∆Φ) , (25)
where ∆x = 2x, ∆Φ = δq • r stands for the relative geometrical phase between the condensate field and its mirror symmetric (originating from the non-zero relative transverse wave-vector δq between them) and I(x) = |E(x, t)| 2 τ for the time-averaged intensity distribution of the sample emission at position x. Here, • • • τ is a time averaging over τ arising from the fact that the camera integration time τ = 1 s is much longer than all time scales involved in the condensate dynamics. In the main text, we implicitly assume that the ergodic hypothesis is valid, which implies that averaging physical observables over long time (as done experimentally) or over a large set of different noise realizations (as in simulations) is equivalent. In what follows, • • • indistinctly denotes temporal or statistical averaging.

The first order correlation function g (1) in Eq. ( 25) is defined by: 27). The photoluminescence background strongly interferes with itself at (∆x = 0, ∆t = 0), creating a peak in |g (1) | at ∆x = 0. This peak disappears after one scan step only, as the photoluminescence background light is fully incoherent.

g (1) (∆x, ∆t) = E (x, t 0 )E(-x, t 0 + ∆t) |E(x, t 0 )| 2 |E(-x, t 0 + ∆t)| 2 . ( 26
)
Experimentally, we retrieve the correlation function ( 26) by measuring the fringe visibility:

V = (I + -I -)/(I + + I -)
, where I + and I -stand respectively for the upper and lower envelopes of I c . At every time delay ∆t, I + and I -are extracted using Fourier analysis on the interferogram. As an example, Fig. 8b and Fig. 8c respectively show the lower and upper envelopes associated to the interference pattern in Fig. 8a. A cut of those three intensity maps along x = 0 yields the graph in Fig. 8d. The visibility V is finally related to the first-order coherence through: (1) | in linear scale. Errorbars are calculated by performing a repeatability analysis on the numerical extraction of g (1) from the interferograms. We report on this plot the result of the fitting procedure describe in III E 2 (red lines).

V (∆x, ∆t) = 2 I(x)I(-x) I(x) + I(-x) |g (1) (∆x, ∆t)| = K(∆x) |g (1) (∆x, ∆t)|, (27)

Estimation of the scaling exponents β and χ

In the main text, we show that the temporal (∆x = 0, (1) | qualitatively agree with the stretched exponential scaling predicted by KPZ theory. In this section, we present a more quantitative analysis of the experimental data, based on curve fitting, which aims at measuring the universal scaling exponent χ and β within a 95% confidence interval.

As mentioned in the main text, the theoretical value of the roughness exponent χ = 1/2 is not characteristic to the KPZ universality class but rather shared among three different classes: Edward-Wilkinson (χ = 1/2, β = 1/4), KPZ (χ = 1/2, β = 1/3) and the class (χ = 1/2, β = 1/2) to which linear systems described by Bogoliubov theory pertain. As all the classes to which our system could belong share the same value for χ, we first set χ = 1/2 and fit the spatial decay of |g (1) (∆x, ∆t = 0) | with the stretched exponential function f x (∆x) = exp (-∆x/λ) /κ (see Fig. 10a). The normalization factor κ and the nonuniversal space-scale λ are two fitting parameters. From this fit, we obtain: κ = 1.14 ± 0.01 (uncertainties are estimated from the 95% confidence interval on the fit parameters). We then focus on the temporal decay of |g (1) (∆x = 0, ∆t) |. In order to properly propagate the error on κ, we first renormalize the data points by defining g (1) = κg (1) . If we omit the uncertainty on the experimental data and only consider the uncertainties originating from the fitting procedure, the error on g (1) can simply be expressed as: δg (1) = δκ g (1) . The renormalized data points are shown on Fig. 10b (on this graph, errorbars are smaller than the points diameter). We finally fit the temporal decay of |g (1) (∆x = 0, ∆t) | with an other stretched exponential function f t (∆t) = exp -(∆t/τ ) 2β , using a weighted nonlinear least squares algorithm to take the error on g (1) into account. From this fitting procedure, we obtain: β = 0.35 ± 0.02. The fitted value of β is in close agreement with the theoretical prediction for the KPZ universality class where β th = 1/3. We can push our analysis further relaxing the constraint on χ. We fit the spatial decay of |g (1) (∆x, ∆t = 0) | with the stretched exponential function fx (∆x

) = exp [-(∆x/λ) 2χ ] /κ,
where χ is now a third fitting parameter (see Fig. 11a). We obtain: κ = 1.14 ± 0.06 and χ = 0.51 ± 0.08. Propagating the error on κ in the same way as before and fitting the renormalized data points by f t (see Fig. 11b), we get: β = 0.36 ± 0.11. As expected, the uncertainty on β is now larger but it still allows us to discriminate between the different universality classes, as the KPZ value β th = 1/3 remains the only one lying within the 95% confidence interval on β. (1) | at ∆t = 0 by the stretched exponential function fx . b Fit of the temporal decay of |g (1) | = κ|g (1) | at ∆x = 0 by the stretched exponential function f t . Errobars are estimated by propagating the uncertainty on κ, as explained in the text. In both graphs, the red dots indicate where the fits are performed. The grey-shaded area gives the 95% confidence interval on those fits.

Symmetric Lieb lattice

As mentioned in section III A, the cavity-exciton detuning of the symmetric lattice is 3 meV larger (in absolute value) than the asymmetric lattice one. As a consequence, the interplay between gain and dissipation gives rise to condensation in the P-bands of the symmetric lattice, at an energy close to the one at which condensation was observed in the asymmetric one. The black arrows in the low-power far-field photoluminescence (see Fig. 12a) indicate the top of the P-band in which the polariton condensate forms. The polariton mass is negative there, thus preventing the formation of modulation instability.

Note the distinctive spatial distribution of the condensate (visible on the interferogram in Fig. 12d), that exhibits two lobes on each pillar, confirming the fact that condensation occurs in P-bands.

Following the procedure described in section III D 2, we can retrieve |g (1) (∆x, y)| for any time delay ∆t and then study the spatio-temporal scaling of -2log |g (1) | for the symmetric lattice. The experimental data obtained for P/P th ≈ 1.12 are shown in Fig. 13. The variations of -2log |g (1) | as a function of ∆x for ∆t = 0 (Fig. 13a) exhibit a linear trend over the spatial window 15 µm < ∆x < 50 µm (grey shaded area), in agreement with KPZ This result is supported by the observation of a plateau in D x = -2∂log(|g (1) (∆x, 0)|)/∂∆x (inset) in the same range of ∆x. The variations of -2log |g (1) | as a function of ∆t 2/3 for ∆x = 0 (Fig. 13b) clearly show a linear increase over the temporal window 25 ps < ∆t < 90 ps (grey area), indicating that this quantity scales as a ∆t 2β power law, with β = 1/3. Finally, we observe in Fig. 13c the collapse onto a single curve of all the data points lying within the non-hatched region of |g (1) | (see inset). This curve can be reproduced with remarkable agreement using the KPZ universal scaling function (black solid line), that has been shifted horizontally and vertically to fit the data points. These results highlight the fact that our experimental findings apply to different lattices and different types of bands.

Effect of the pumping power on the KPZ window

We briefly discuss in this section the impact of the excitation power on the coherence decay of the polariton condensate. Fig. 14 shows, for ∆x = 0, the variations of -2log(|g (1) |)

as a function of ∆t 2/3 for different excitation powers, both on the asymmetric (left) and symmetric (right) Lieb lattice. The linear trend observed at low excitation power (P/P th < 1.15), emphasized by the blue lines on both panels, becomes less and less visible when P/P th increases. Moreover, the extension of the temporal KPZ window over which this linear trend occurs shrinks progressively as the power increases. Simulations presented in Fig. 14c reproduces the observed features and are discussed in Sec. V B 4. The power dependent behavior is further supported by the results shown in Fig. 15, where the collapse of the data onto the KPZ scaling function is plotted for three values of the power (P/P th = 1.13, 1.20 and 1.26) used to illuminate the asymmetric Lieb lattice. The spatio-temporal KPZ window (non-hatched region in Fig. 15a2, Fig. 15b2 and Fig. 15c2), defined by the data points in the |g (1) | data set that collapse onto F KPZ (black solid line), shrinks as we increase the power. reproduce the spatial density profile of the condensate in our experiments, we use a spatiallydependent pump profile modelled by

P (x) = P [1 + tanh((L 0 + x)/σ)] [1 + tanh((L 0 -x)/σ)] [1 + tanh(L 0 /σ)] 2 (31)
with L 0 = 80 µm and σ = 9.7 µm. The remaining free parameters in our numerical simulations are thus the scattering rate R of excitons into the condensate and the reservoir decay rate γ R . We adjust them so as to find the best possible agreement with the experimental data.

All the simulations presented in this article were performed using R = 8.8 × 10 -4 µm ps -1

and γ R = 0.45γ 0 .

In connection with the derivation of the mapping from the gGPE equation for the condensate field to the effective KPZ equation for the phase field presented in Sec. II B, let us comment on the involved time scales. It is not straightforward to determine the dynamical time scales for density and phase fluctuations which are in general due to many-body non-linear effects. However, one can assume that they are typically set by the one-body relaxation rates γ R and γ 0 for the densities, and by the mean field frequency ω 0 = 2g R n R + gρ 0 for the phase. Thus one has ω 0 12γ 0 6γ R , which suggests that the decoupling of time scales assumed in the mapping to the KPZ equation is verified. In our system, the phase dynamics is faster than the density dynamics of both the condensate and the reservoir.

Moreover, the numerical simulations confirm the emergence of KPZ dynamics for our set of parameters.

This mapping allows us to evaluate the values of the parameters of the KPZ equation.

Using the previous values for the gGPE parameters, we find:

λ = -5.7 × 10 2 µm 2 ps -1 , ν = 3.8 × 10 2 µm 2 ps -1 , D = 2 µm . (32)
This allows us in the following to compute the theoretical values for the non-universal normalization constants entering the universal scaling function and distribution, and compare them with the results from direct fits of the experimental and numerical data.

Note that in the main text, the black curve that we adjust to the data points is the KPZ universal scaling function for the stationary case, which is the only one that is known exactly [START_REF] Prähofer | Exact scaling functions for one-dimensional stationary kpz growth[END_REF]. FIG. 16. Collapse of the simulated data. Collapse of the numerical data for |g (1) (∆x, ∆t)| plotted in dimensionless units, using the normalization constants y 0 and Γ determined from the microscopic coefficients of Eq. [START_REF] Bobrovska | Stability and spatial coherence of nonresonantly pumped exciton-polariton condensates[END_REF]. The dashed lines show the exact values of horizontal asymptotes expected for the flat, curved and stationary universality subclasses.

B. Different contributions to g (1) and influence of space-time vortices

In the previous sections, we have demonstrated that both the experimental and the numerical data for the first-order coherence g (1) exhibit the KPZ scaling in space and time, yielding a clear collapse of all data points within the KPZ window onto the universal scaling function. The crucial question which arises is to what extent the behavior of g (1) properly reflects the properties of the phase itself. This is all the more important that, as shown in the main text, typical space-time phase maps exhibit the formation of space-time vortices.

In this section, we analyze the effect of these vortices and show that (i) they do not spoil the KPZ regime as long as almost all vortices appear as pairs of close-by vortex and anti-vortex and the density of single vortices remains low enough, and (ii) the scaling behavior of g (1) is indeed inherited from the scaling behavior of the phase-phase correlations.

Analysis of the effect of density-density and density-phase correlations

We first examine the different contributions to g (1) in order to test the assumptions made in Sec. II C to relate the first-order coherence of the condensate field to the phase-phase correlations. The assumptions required to derive Eq. ( 19) from Eq. ( 15) are to neglect both the density-phase and density-density correlations. The calculated temporal variations of these correlations are displayed in Figs. 17a andb, for ∆x = 0. We observe that although such correlations are present in the system, they remain approximately constant over the KPZ window. Furthermore, Fig. 17c shows the effect of each of these contributions, by comparing:

(1) The quantity -2 log (| exp(i∆θ) |) defined in Eq. ( 19) when only considering phasephase correlations in g (1) (black line),

(2) The right hand side of Eq. ( 16), obtained by considering phase-phase and densitydensity correlations in g (1) (blue line), 16), neglecting density-phase correlations in g (1) . The black line shows the variations of -2 log(| exp(i∆θ) |) where density-density and density phase correlations are neglected. The grey shaded region in all panels indicates the KPZ window.

we are able to find a vortex-free region that extends over a time interval exceeding 100 ps.

We thus use these trajectories without any further processing for evaluating the vortex-free phase variance. For about 20 % of the realizations, any 100 ps time window contains at least one or a few V-AV pairs, but that are sparse enough to allow the numerical filtering of the phase jumps and so to evaluate the vortex-free phase variance (in practice, we filter a vortex by adding at its location a second vortex of opposite charge). In a limited number of cases (5 % of the trajectories) the jumps are too numerous in any time window of 100 ps to allow computing the vortex-free variance. We thus discard these trajectories in the calculation of the phase variance.

The computed vortex-free phase variance is shown in Fig. 17c (yellow diamonds) and is directly compared with -2 log e i∆θ . Both quantities perfectly coincide within the KPZ window, up to the departure from this regime, which shows that the cumulant expansion Eq. ( 21) is valid in this whole range. We emphasize that -2 log e i∆θ is computed over the full duration of all unprocessed trajectories, and thus includes the effect of vortices.

Therefore, we conclude that in the low-vortex density regime we explore here, the quantity -2 log e i∆θ (and thus also g (1) ) is not sensitive to V-AV pairs. This analysis fully confirms that g (1) provides a good observable to probe the KPZ scaling of the phase. residual linear behavior becomes dominant compared to the ∆t 2/3 scaling of the fluctuations and introduces faster decay of the coherence. This demonstrates that the departure from the KPZ regime in the raw correlations is induced by inhomogeneous spectral broadening.

Note that we also performed simulations for higher power densities. As reported in Fig. 14c, we observe that departure from the KPZ regime occurs at shorter time delays when increasing the excitation power, thus suggesting that inhomogeneous broadening increases with excitation power. We point out that the shrinking of the KPZ window obtained in the numerical simulations occurs over a range of excitation powers that is similar to the one observed experimentally.

C. Distribution of phase fluctuations

In this section, we discuss the probability distribution associated to phase fluctuations computed from our numerical simulations. We demonstrate that this distribution is well reproduced by the Tracy-Widom (TW) Gaussian Orthogonal Ensemble (GOE) distribution, which is characteristic of the flat KPZ universality subclass.

Let us emphasize that by normalizing the variance of the distribution at each ∆t to unity, one can deduce the estimate for the non-universal parameter |Γ| as

|Γ num | = 1 ∆t
Var(∆θ(∆t))

Var T W GOE 3/2 . ( 38
)
The values found for the different time instants can be compared with the theoretical estimate of |Γ| given in Eq. [START_REF] Fischer | Spatial coherence properties of one dimensional exciton-polariton condensates[END_REF]. The comparison is shown in the right panel of Fig. 20, and the agreement is remarkable. This is yet another strong confirmation of the accuracy of the mapping, given that Γ is a non-universal constant. [START_REF] Fischer | Spatial coherence properties of one dimensional exciton-polariton condensates[END_REF]. It is important to note that numerically we only have access to the absolute value of Γ because it is extracted from the variance of the phase, which is a positive defined quantity.

  sample consists of a high quality factor (Q 70000) λ/2 Ga 0.05 Al 0.95 As microcavity surrounded by two Al 0.20 Ga 0.80 As/ Al 0.05 Ga 0.95 As distributed Bragg reflectors. Three stacks of four 8 nm GaAs quantum wells are embedded in this microstructure, at the anti-nodes of the cavity mode electromagnetic field, resulting in a 15 meV collective Rabi splitting. The as-grown planar cavity is patterned into 1D lattices of coupled micropillars (3 µm in diameter), using electron beam lithography and dry etching. In this work, we use a 200 µm long asymmetric Lieb lattice, made of three pillars per unit cell with 2.2 µm center-to-center separation distance. A close-cycle cryostation holds the sample at 4K.

FIG. 1 .FIG. 2 .

 12 FIG. 1. KPZ physics in the phase dynamics of a 1D polariton condensate. a Snapshots of the condensate phase evolution obtained by numerically solving Eqs. (2)-(3). The unwrapped phase evolves in time as a growing interface. b Intensity distribution of the condensate emission measured for P/P th = 1.1. c Sketch of the lattice together with the excitation scheme. The lattice is excited using an elongated flat top beam. Inset: Sketch of a micropillar inner structure. A semiconductor optical microcavity is enclosed by two Bragg mirrors. Quantum wells (QWs, red layers) are distributed at the anti-nodes of the cavity mode (black line).

10

 10

FIG. 3 .

 3 FIG. 3.KPZ scaling in the coherence decay of a 1D polariton condensate. a Measured values of -2 log |g(1) (0, ∆t)| as a function of ∆t 2/3 . Inset: Temporal derivative of -2 log |g(1) (0, ∆t)| , computed from the experimental data, as a function of ∆t. b Measured values of -2 log |g(1) (∆x, 0)| as a function of ∆x. Inset: Spatial derivative of -2 log |g(1) (∆x, 0)| , computed from the experimental data, as a function of ∆x. On both panels, the grey area delimits the KPZ temporal window. In panels (a) (resp. (b)), the simulated data are shown with an orange line (resp. squares). Errorbars on the experimental data points are calculated by performing a repeatability analysis on the numerical extraction of g(1) (∆x, ∆t) from the interferograms. The orange-shaded area (resp. orange errorbars) on panel (a) (resp. (b)) show the 95% confidence interval on the simulated data. c Measured values of -2 log(κ|g(1) |)/∆t 2/3 as a function of y = ∆x/∆t 2/3 , for points within the non-hatched region of the coherence map (bottom inset). Top inset: Same as in the main panel but for numerical data. On both graphs, the black line corresponds to the KPZ scaling function F = C 0 F KPZ (y/y 0 ), adjusted to the experimental data by tuning the values of C 0 and y 0 . Experimental and simulated |g(1) | datasets are normalized using κ = 1.13. In all panels, the numerical data are averaged over 10 4 realizations of the noise.

FIG. 4 .

 4 FIG. 4. Analysis of the simulated phase dynamics. a Example of a phase map obtained after subtraction of the dynamical phase, ω 0 ∆t. Inset: Zoom on a V-AV pair. b Time evolution of the unwrapped phase along x = 0 (blue line): two V-AV pairs are crossed, resulting in two phase jumps of approximately 2π. The black line of panel (b) is obtained by filtering all the vortices present in the phase map (a). c -2 log |g (1) | (red line) and Var [∆θ VF ] (symbols) as a function of ∆t 2/3 for ∆x = 0, computed over 10 4 realizations. The orange line is the same as in Fig. 3a. Inset: Temporal derivative of -2 log |g (1) | (red line) and Var [∆θ VF ] (symbols) as a function of ∆t, showing similar behaviour in the KPZ window. d Probability distribution of phase fluctuations obtained at ∆t = 50 ps (blue squares) from numerical simulations. Our computation involves 3 × 10 4 realizations of the noise. The first (resp. second, third, etc) peak from the left gathers the realizations for which no jump (resp one, two, etc) has occurred before ∆t = 50 ps. Yellow line: Fit of the numerical data by the Tracy-Widom GOE distribution.

F

  .B. contributed to the design of the sample structure and initial characterization of the sample. A.L., M.M. grew the sample by molecular beam epitaxy. I.S, L.L.G., and A.H. fabricated the polariton lattices. Q.F., D.S., I.A., M.W., I.C., A.A., M.R., A.M., L.C., S.R., and J.B. participated to the scientific discussions about all aspects of the work. Q.F., A.M., L.C., S.R. and J.B. wrote the original draft of the paper. Q.F., D.S., I.A., M.W., I.C., A.A., M.R., A.M., L.C., S.R., and J.B. reviewed and edited the paper into its current form. A.M., L.C., S.R. and J.B. supervised the work.

FIG. 5 .

 5 FIG. 5. Sketches of the lattices. a Asymmetric and b symmetric Lieb lattices.

FIG. 6 .

 6 FIG. 6. a Dispersion of the upper S-band obtained by fitting the spectrum (inset) by a lorenztian lineshape for all wave-vectors lying in the first Brillouin zone. Red dotted line: parabolic fit of the data points in the vicinity of k = 0. Yellow line: fit of the dispersion using a nonlinear regression model. b Group velocity as a function of k, computed from the data in (a). Yellow line: derivative of the fit function in (a). c Spectral linewidth γ spe as a function of k.

FIG. 7 .

 7 FIG. 7. a Energy-resolved photoluminescence in real space (log scale), obtained under local nonresonant excitation in the vicinity of x = 0. Bottom inset: Exponential decay of the luminescence intensity at E = 1.5673 meV (log scale). Fitting this decay (black solid line) enables to retrieve the Lorentzian contribution to γ spe . b Lorentzian (γ pro , red diamonds) and Gaussian contribution (γ g , blue dots) to γ spe . c Spectral linewidth γ spe (light blue dots) and polariton linewidth γ (red stars) as a function of k. The red dots are obtained from the blue ones by removing the inhomogeneous Gaussian broadening in γ spe . d Polariton linewidth γ together with our theoretical model (blue line). The dashed red line shows the quadratic expansion of the model in the vicinity of k = 0.

1 .

 1 Optical setupThe sketch of the optical setup is shown on Fig.2dof the main text. In our experiment, polaritons are excited using a non-resonant continuous-wave laser of wavelength 740 nm (where the cavity mirror reflectivity exhibits a minimum). A spatial light modulator (SLM) enables shaping the excitation spot into a 125 µm long flat-top beam in the lattice direction, and a Gaussian with a 3.5 µm FWHM in the transverse direction. The light emitted by the sample is collimated by the excitation lens, passes through a polarizer (selecting the TM polarization), and is sent through an interferometer. A polarized beam splitter, combined with a half wave-plate, enables splitting the incoming light into two beams while controlling their power ratio. The first beam reflects on a plane mirror, making a round trip through a quarter-wave plate which turns its polarization by 90 o . The second beam reflects on a retroreflector. The latter is mounted on a motorized translation stage, allowing for a variation of the path length difference between the interferometer arms. Both beams are finally recombined in a non-polarized beam splitter before being imaged onto a CCD camera.

FIG. 8 .

 8 FIG. 8. Retrieving the first order coherence of polariton condensates. a Interference pattern captured by the CCD for ∆t = 0. b-c lower and upper envelopes of the interferogram. d Cut of the first three maps along x = 0. Coherence retrieved from Eq. (27). The photoluminescence background strongly interferes with itself at (∆x = 0, ∆t = 0), creating a peak in |g(1) | at ∆x = 0. This peak disappears after one scan step only, as the photoluminescence background light is fully incoherent.

1 FIG. 9 .

 19 FIG.9. Variations of |g(1) | in linear scale. Errorbars are calculated by performing a repeatability analysis on the numerical extraction of g(1) from the interferograms. We report on this plot the result of the fitting procedure describe in III E 2 (red lines).

  Fig 3.a) and spatial variations (∆t = 0, Fig 3.b) of |g

FIG. 10 .

 10 FIG.10. Fit of the coherence spatial and temporal decay setting χ = 1/2. a Fit of the spatial decay of |g(1) | at ∆t = 0 by the stretched exponential function f x . b Fit of the temporal decay of |g(1) | = κ|g(1) | at ∆x = 0 by the stretched exponential function f t . Errobars are estimated by propagating the uncertainty on κ, as explained in the text. They are smaller than the points diameter. In both graphs, the red dots indicate where the fits are performed. The grey-shaded area gives the 95% confidence interval on those fits. Inset: Residual Sum of Squares (RSS) as a function of β, reaching a minimum at 0.35.

1 FIG. 11 .

 111 FIG.11. Fit of the coherence spatial and temporal decay setting χ as a fitting parameter. a Fit of the spatial decay of |g(1) | at ∆t = 0 by the stretched exponential function fx . b Fit of the temporal decay of |g(1) | = κ|g(1) | at ∆x = 0 by the stretched exponential function f t . Errobars are estimated by propagating the uncertainty on κ, as explained in the text. In both graphs, the red dots indicate where the fits are performed. The grey-shaded area gives the 95% confidence interval on those fits.

FIG. 12 .

 12 FIG. 12. Condensation in the symmetric Lieb lattice. a-b Far-field energy-resolved photoluminescence (TM polarization) below (a) and above (b) condensation threshold. The condensate forms at the edge of the second Brillouin zone, in a negative mass state lying in the P-bands of the Lieb lattice. c Integrated emission intensity as function of the incident pump power, showing a condensation threshold at P th ≈ 28 mW. d Interference pattern at ∆t = 0. The black circles show where the pillars are located.

4 FIG. 13 .

 413 FIG. 13. KPZ scaling in the coherence decay of a polariton condensate forming in the symmetric Lieb lattice. a Measured values of -2 log(|g (1) (∆x, 0)|) as a function of ∆x. The grey shaded area highlights the temporal window where KPZ scaling settles in. Inset: Variations of D x computed from the experimental data. The solid line shows the expected scaling of D x . b Measured values of -2 log(|g(1) (0, ∆t)|) as a function of ∆t 2/3 . The grey shaded area shows the temporal window inside which KPZ scaling is observed. Errorbars on the experimental data points are calculated by performing a repeatability analysis on the numerical extraction of g(1) (∆x, ∆t) from the interferograms. c Measured values of -2 log(κ|g(1) |) as a function of the rescaled coordinates y = ∆x/∆t 2/3 . We observe the collapse of the data points within the non hatched region of the coherence map (inset) onto the KPZ universal scaling function (black solid line). The normalization factor is κ = 1.32.

FIG. 14 .

 14 FIG.14. a-b Measured values of -2log(|g(1) (0, ∆t)|) as a function of ∆t 2/3 for various excitation powers in (a) the asymmetric and (b) the symmetric Lieb lattice. For P/P th < 1.15, the data exhibits a linear increase over a given time window characteristic of KPZ scaling. This window shortens as the excitation power increases before vanishing completely. c Computed values of -2log(|g(1) (0, ∆t)|) as a function of ∆t 2/3 , for various excitation powers (p = P/P th varies from 1.13 to 1.25, as indicated on the curves). Our simulations qualitatively reproduce the change observed experimentally in the scaling of -2log(|g(1) (0, ∆t)|) when increasing the excitation power. The parameters used in our numerical analysis are detailed in section IV.

FIG. 15 .

 15 FIG. 15. Shrinking of the spatio-temporal KPZ window when increasing the excitation power. a P/P th = 1.13. b P/P th = 1.20. c P/P th = 1.26. The normalization factor κ is equal to 1.13 in (a-b) and to 1.11 in (c).

- 2 FIG. 17 .

 217 FIG.17. a Density-density correlations. b Density-phase correlations. c Successive approximations of -2 log(|g(1) |). The red line is the full reconstruction of -2 log(|g(1) |), taking into account phase-phase, density-density and density-phase correlations. The blue one is obtained from Eq. (16), neglecting density-phase correlations in g(1) . The black line shows the variations of -2 log(| exp(i∆θ) |) where density-density and density phase correlations are neglected. The grey shaded region in all panels indicates the KPZ window.

FIG. 18

 18 FIG. 18. a, b and c. Set of three space-time phase maps generated in the numerical simulations using the parameters given in section IV, with (a) no vortices, (b) few vortices, (c) many vortices. The most likely realization is the one depicted in panel (b). d Corresponding cuts at x = 0, showing that the phase trajectory at x = 0 undergoes jumps whenever it crosses a vortex-anti-vortex pair.

5 FIG. 19 .

 519 FIG. 19. a Set of unwound phase trajectories with no phase jump, after subtracting the ensemble average linear growth ω 0 ∆t. Inset: Distribution of the mean frequency ω 0 of each individual trajectory. The fluctuations of the frequency ω 0 around its ensemble average value ω0 induces a residual slope in these trajectories. b Correlations -2 log | e i∆θ | with and without correcting for the frequency fluctuations.

  FIG.[START_REF] Krug | Origins of scale invariance in growth processes[END_REF]. a Distribution of the phase fluctuations obtained by summing up the individual distribution obtained at all delays within the range 30 ps ≤ ∆t ≤ 80 ps. The Tracy Widom GOE distribution is shown as a yellow solid line. For comparison, we also plot with green solid line the normal distribution N , centered at ∆θ = 0 with the same variance as the Tracy-Widom distribution. In the inset, the agreement between the numerical data and either (green) the Normal distribution or (yellow) the Tracy-Widom distribution are compared by showing the ratio P[∆θ]/P α [∆θ] with α = T W GOE , N . This clearly demonstrates that TW GOE provides a better description of the numerical data. b Distribution of the parameters |Γ num | in the same temporal window, extracted from the simulations via (38), together with the theoretical value |Γ th | computed using the microscopic parameters of the numerical simulation and[START_REF] Fischer | Spatial coherence properties of one dimensional exciton-polariton condensates[END_REF]. It is important to note that numerically we only have access to the absolute value of Γ because it is extracted from the variance of the phase, which is a positive defined quantity.

Acknowledgments. We thank Valentin Goblot, Daniel Vajner and Ateeb Toor for their assistance in the early development of the experiment.

Funding. This work was supported by the Paris Ile-de-France Région in the framework of DIM SIRTEQ, the French RENATECH network, the H2020-FETFLAG project PhoQus (820392), the QUANTERA project Interpol (ANR-QUAN-0003-05), the European Research

Supplementary Material

where K(∆x) is a normalization factor taking into account potential imbalance between I(x) and I(-x). In our case, this factor remains close to 1. Using Eq. ( 27), we finally retrieve |g (1) (∆x, y; ∆t)| (see Fig. 8e). The coherence map shown in Fig. 2e of the main text is obtained by keeping only the maximum value of |g (1) (∆x, y; ∆t)| over the pillars identified through white solid circles in Fig. 8e. (1) | After having retrieved |g (1) (∆x, ∆t)| from the data analysis detailed in the previous section, we search for KPZ scalings in the spatio-temporal variations of -2log |g (1) | . In particular, we show in Fig. 3c of the main text the collapse onto the universal KPZ scaling function of the |g (1) | data points within a certain spatio-temporal window. In order to do so, we plot in log-log scale -2log κ|g (1) | /∆t 2/3 as function of the rescaled coordinate y = ∆x/∆t 2/3 , where κ is a normalization factor that needs to be properly set. Indeed, representing the data in such a way implicitly requires that the temporal KPZ scaling extends all the way to ∆t = 0, where |g (1) (0, 0)| is expected to be 1. Experimentally, we observe a transient regime at short time delays, preceding the establishment of the ∆t 2β power law behavior of -2log |g (1) | . We thus need to ensure that the extrapolation of this power law passes through 0 at ∆t = 0 in order for the chosen graphic representation to be meaningful. This amounts to shifting downward the data points shown in Fig. 3 a-b until they match the blue solid line, which, in turn, translates into multiplying the whole |g (1) | data set by a factor κ. Note that this normalization does not change the coherence decay, which remains a stretched exponential. (1) | in linear scale

Normalization of |g

E. Additional results

Variations of |g

In order for the reader to be readily able to compare our results with with other works in the literature where data are reported in a different way, we show in figure 9 the variations of |g (1) | as a function of ∆t (for ∆x = 0) and ∆x (for ∆t = 0). We also report on this graph the fit of the data points by streched-exponential decays over the spatial and temporal KPZ window. Further details on the fitting procedure can be found in section III E 2.

IV. NUMERICAL SIMULATIONS: METHOD AND PARAMETERS

The numerical integration of Eq. ( 6) is performed using the interaction picture method [START_REF] Werner | Robust algorithms for solving stochastic partial differential equations[END_REF][START_REF] Dennis | Xmds2: Fast, scalable simulation of coupled stochastic partial differential equations[END_REF]. The idea behind this integration scheme is similar to the interaction picture in quantum mechanics. We first split Eq. ( 6) into a linear, exactly solvable part and a remaining nonlinear part. We then solve the linear component in Fourier space and transform it back to real space. We transform Eq. ( 6) by moving into the interaction picture and integrate the resulting nonlinear equation using semi-implicit Runge-Kutta method, with an adaptive time-step. We take as initial condition ψ(x, t = 0) = 0, and let the condensate grow under the action of the pumped reservoir. The sampling starts at t 0 , long after the condensate has reached its stationary density profile. Usually, we perform our simulations using t 0 = 10 ns.

As mentioned in section III C, some of the parameters entering the numerical simulations are known experimentally. For instance, we use a lattice spacing equal to the experimental lattice period a = 4.4 µm. The measurement of the polariton dispersion relation E(k) shown in Fig. 6a provides us with a good estimate of the polariton mass m = -3.3×10 -6 m e . The k-dependent polariton linewidth γ(k) is modelled by the function

which is compared to the experimental data in the inset of Fig. 7d (blue solid line). The

. Furthermore, the energy blue-shift at threshold is known with good accuracy:

where n R | P =P th = P th /γ R = γ 0 /R stands for the reservoir density at condensation threshold and P th = γ 0 γ R /R for the threshold power. Assuming that g R does not depend on the pumping power P , we find:

Since the reservoir-induced blueshift 2g R n R is two orders of magnitude larger than the polariton-induced blueshift g|ψ| 2 , we take g = 0 in our simulations. In order to qualitatively

V. NUMERICAL SIMULATIONS: DISCUSSION

A. Deducing the universality subclass from the collapse of numerical data

In this section, we explain how the horizontal asymptote of the curve onto which the numerical data points collapse provide information about the KPZ universality subclass the system belongs to. The KPZ universal scaling function associated with the correlation function Var [∆θ(∆x, ∆t)] in 1D is defined by:

where y 0 and Γ are non-universal normalization constants which can be directly expressed in terms of the KPZ parameters λ, ν and D as follows

Using the parameters given in section IV, we obtain:

In Fig. 16, we report the same data points as shown in the inset of Fig. 3c but using dimensionless coordinates, that we calculate based on the values of Γ and y 0 given above. We recognize the typical features of the KPZ scaling function, showing a plateau at small y and a linear growth at large y [START_REF] Deligiannis | Accessing kardar-parisi-zhang universality sub-classes with exciton polaritons[END_REF]. It is worth mentioning that the horizontal asymptote F KPZ (0) of the universal scaling function is a universal constant that depends on the universality subclass the system falls into. It is known exactly for the flat, stationary and curved initial conditions (see Sec. V C for details). In Fig. 16, we indicate by horizontal dashed lines the values of F KPZ (0) expected in the three different KPZ subclasses. We observe that the plateau reached by the simulated data at small y is in close agreement with the exact theoretical value for the flat subclass F KPZ,flat (0) 0.64 (which is the one expected in our case, as we will see Sec. V C). This agreement is remarkable given that Γ and y 0 are non-universal quantities (i.e. sensitive to the microscopic details of the model). It provides even at a quantitative level a strong support for the validity of our mapping between the gGPE model and the KPZ equation for the phase dynamics (see Sec. II B).

Comparing the black and blue curves, one observes that except at very short time delays (∆t < 3 ps), density-density correlations only lead to a global shift, and thus do not affect the scaling behavior. When including density-phase correlations (red curve) the main effect we observe is a faster decoherence at long time delays (∆t > 80 ps). Overall, Fig. 17 shows that the scaling of the different computed quantities stays unaffected within the KPZ temporal window. From this analysis, we conclude that Eq. ( 19) is a reliable approximation of g (1) in the KPZ regime, which amounts to saying that the behavior of g (1) is mainly dominated by phase-phase correlations in the KPZ window. In the following, we finally relate e i∆θ to the variance of the phase Var [∆θ].

Spatio-temporal phase maps and calculation of the phase variance

A visual inspection of typical space-time phase maps reveals the presence of space-time vortices. Fig. 18 shows three such maps corresponding to three different noise realizations presenting no (a), few (b) or many (c) space-time vortices. By analyzing the distribution of typical vortex distances, we have found that almost all space-time vortices appear as vortex-antivortex (V-AV) pairs (set of two nearby vortices of opposite charge), while the number of single vortices is negligible. The associated cuts at x = 0 are shown in Fig. 18d.

When no vortices are present, the unwrapped phase at x = 0 shows a dominant linear behavior in time (dark blue line in Fig. 18d), on top of which KPZ fluctuations develop.

When crossing V-AV pairs, the unwrapped phase at x = 0 undergoes jumps on very short timescales (light blue line in Fig. 18d). Every jump induces a phase shift with respect to vortex-free trajectories, after which the linear behavior is restored. The amplitude of these jumps is distributed between 0 and 2π depending on where the x = 0 line crosses the V-AV pair, but dominated by values close to 2π.

Clearly, such phase jumps will have a strong impact on the calculation of the variance of phase fluctuations: vortices will lead to a fast increase of the variance that may hide the KPZ scaling, or even lead to other dynamical regimes such as the ones evidences in Ref. [START_REF] He | Space-time vortex driven crossover and vortex turbulence phase transition in one-dimensional driven open condensates[END_REF]. For most of the phase maps generated in the simulations with realistic experimental parameters, we notice that only few jumps occur within the 300 ps time interval under consideration.

More quantitatively, over a set of 10 4 trajectories, the probability of observing a jump within a 1 ps time interval is around 0.01. As a consequence, for most realizations (about 75 %),

Resilience of KPZ to space-time V-AV pairs

In the previous Section, we have found that the presence of V-AV pairs weakly affects the temporal scaling of |g (1) (0, ∆t)|, as can be observed in Fig. 17. The analysis of the phase jump amplitudes in Fig. 4b of the main text provides an explanation for the robustness of the g (1) correlations against space-time V-AV pairs. Indeed, we notice that the jumps are centered around values that are close to multiples of 2π. Hence, such jumps have a negligible impact as they enter in the exponential of the condensate field. This property allows observing the emergence of KPZ universality even in the presence of some space-time V-AV pairs. This property, characteristic of the compact version of the KPZ problem, is even more remarkable if we notice that, in such conditions, Eq. ( 21) cannot be blindly used:

in presence of defects every random jump increases considerably the phase-phase correlator bringing it away from the predictions of the KPZ scaling.

Effect of the condensate linewidth on the long-time coherence

For exciton-polaritons, the mean velocity ω 0 of each phase trajectory fluctuates with respect to the ensemble average ω 0 . To illustrate this, we show in Fig. 19 a set of trajectories with no jumps in the KPZ time window, where we subtracted the average dynamical phase ∆θ(0, ∆t) ≡ ω 0 ∆t. One observes a dispersion of the trajectories, which are not on average constant in time but rather exhibit a residual linear behavior, with slope δω 0 = ω 0ω 0 . This is because the total frequency ω 0 of each individual trajectory varies due to density fluctuations. Such stochastic variations in the slope of the phase dynamics corresponds to an inhomogeneous spectral broadening. The related distribution P(ω 0 ) is shown in the inset of Fig. 19a.

To determine the influence of this broadening on the KPZ regime, we post-process the phase trajectories so as to subtract for each trajectory its own linear behavior ω 0 ∆t. We then compute the corresponding correlations | e i∆θ-ω 0 ∆t |. The result is shown in the right panel of Fig. 19, where it is compared with the "raw" | e i∆θ | where the average is performed over the same set of trajectories, but subtracting the average dynamical phase ω 0 ∆t. This clearly evidences that the KPZ ∆t 2/3 scaling holds over much longer times once the intrinsic slope of each trajectory is properly removed. Otherwise, at large times (∆t > 80 ps), the

The KPZ universality subclasses

While the KPZ universality class is fully characterized by the critical exponents χ and β for the two-point correlators, the distribution of the height fluctuations of a KPZ interface in 1D allows one to distinguish three universality subclasses. Let us briefly review this result in this section, before detailing our results for the distribution of the phase fluctuations of the exciton-polariton condensates.

For a classical interface, the height field h(t) ≡ h(x 0 , t) is known to behave at long times and at any given point x 0 in space according to

where, v ∞ represents the mean velocity of the growing interface [57], and h is a random variable describing the reduced fluctuations whose distribution is universal. Much theoretical effort has been dedicated during the last decade to study the properties of the dimensionless random field h, whose distribution P [ h] turns out to be sensitive to the spatial profile of the initial condition h(x, t = 0), or equivalently to the global geometry of the interface (we refer to Ref. [START_REF] Takeuchi | An appetizer to modern developments on the kardar-parisi-zhang universality class[END_REF]) for a review). Three main possible subclasses emerge. For flat initial conditions the reduced KPZ field h is distributed according to the Tracy-Widom distribution (TW) associated to the largest eigenvalue of random matrices belonging to the Gaussian Orthogonal Ensemble (GOE). For curved initial conditions, P [ h] is the TW distribution associated to the largest eigenvalue of random matrices belonging to the Gaussian Unitary Ensemble (GUE). The last subclass is associated to stationary initial condition for the phase profile, for which P [ h] is the Baik-Rains distribution.

Based on a previous numerical study [START_REF] Squizzato | Kardar-parisi-zhang universality in the phase distributions of one-dimensional exciton-polaritons[END_REF], we expect that driven-dissipative condensates, in absence of external confinement, belong to the flat universality subclass. This is supported by the universal plateau value obtained in Sec.V A for the KPZ scaling function, which coincides with the one of the flat subclass.

Analysis of the phase distribution

To treat the condensate phase field in a similar way as classical interface height h (unbounded parameter), one should first unwrap the phase in time at fixed x = 0. In the KPZ regime, the equivalent of Eq. ( 36) then becomes:

with ω 0 the mean frequency associated with the phase dynamics and t 0 the reference time for the phase unwinding. We emphasize that the unwrapping is crucial to study the distribution of the reduced variable θ. The distribution of the phase itself has a compact support and cannot be any of the Tracy-Widom or Baik-Rain distributions characteristic of the KPZ realm. Following the analysis usually performed on a classical interface, a natural strategy to obtain the reduced random variable θ from the unwound phase would be to subtract the mean linear behavior and to rescale by ∆t 2/3 . However, as explained in the previous section, the presence of random phase jumps, induced by the proximity of vortices, hinders this strategy, since it disrupts the linear ∆t evolution. In particular, as visible in Fig. 5 of the main text, these phase jumps induce replicas of the main distribution separated by shifts of value close to 2π. In principle, since the amplitude of the phase jumps are distributed, they could also affect the main central peak. However, as shown below, whereas they do affect the right tail of the central distribution, their effect is negligible on its left tail.

To overcome this difficulty, instead of following the time evolution of the unwrapped phase, as explained in the main text, we rather record the distributions of the phase fluctuations over different realizations at fixed time delays ∆t lying in the KPZ window. We then normalize each obtained distribution such that it has zero mean and unit variance. To test with high accuracy the agreement between the phase distribution and Tracy-Widom distribution (in particular to better resolve the tails), we accumulate more statistics by summing the properly normalized fluctuations at all time instants in the KPZ window. The resulting distribution around ∆θ = 0 is displayed in Fig. 20 (left panel). Strikingly, the left tail of the phase distribution strongly departs from a Gaussian distribution and is reproduced using the Tracy-Widom GOE distribution over more than five decades. Note that the right tail cannot be analysed with such precision because it is affected by the subset of trajectories that contain one or more phase jumps.