category: literaturenote citekey: novokreschenovkardarparisizhanguniversality2025 title: Kardar--Parisi--Zhang universality in optically induced lattices of exciton--polariton condensates authors: "Novokreschenov, D.; Neplokh, V.; Misko, M.; Starkova, N.; Cookson, T.; Kudlis, A.; Nalitov, A.; Shelykh, I. A.; Kavokin, A. V.; Lagoudakis, P." year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.10187v1" zotero_key: 7PPSFDZX zotero_storage: JWT7BZ73 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Novokreschenov, D."

status: converted

Kardar–Parisi–Zhang universality in optically induced lattices of exciton–polariton condensates

D. Novokreschenov,¹ V. Neplokh,1, 2 M. Misko,³ N. Starkova,³ T. Cookson,³ A. Kudlis,4, ∗ A. Nalitov,¹ I. A. Shelykh,⁴ A. V. Kavokin,2, 1 and P. Lagoudakis³

(Dated: June 13, 2025)

We investigate space–time coherence in one-dimensional lattices of exciton–polariton condensates formed by fully reconfigurable non-resonant optical pumping. Starting from an open-dissipative Gross-Pitaevskii equation with deterministic reservoir kinetics and stochastic condensate noise, we derive a discrete complex-field model that incorporates coherent tunnelling, reservoir-mediated dissipative coupling and gain-saturation non-linearity. Adiabatic elimination of fast density fluctuations reveals a wedge-shaped region in the complex hopping plane where the coarse-grained phase dynamics reduces to the Kardar–Parisi–Zhang (KPZ) equation. By computing high-resolution phase diagrams of the temporal and spatial scaling exponents we pinpoint the boundaries separating the KPZ domain from the Edwards–Wilkinson (EW) regime. Large-scale graphics processing unit (GPU) simulations of chains containing up to N = 2000 condensates confirm these predictions: inside the wedge the exponents converge to β^N = 0.329(3) ≈ 1/3 and χ^N = 0.504(4) ≈ 1/2, whereas outside it the dynamics moves away from KPZ and ultimately flows toward the EW fixed point, although finite system size and finite observation time may yield intermediate effective exponents. These results pave the way to the implementation of ultrafast KPZ-simulators based on one-dimensional arrays of exciton-polariton condensates.

I. INTRODUCTION

Exciton-polaritons, composite bosonic quasiparticles formed by strong coupling between microcavity photons and quantum-well excitons, constitute a versatile platform for exploring non-equilibrium many-body physics [1]. One fascinating aspect is that polaritons can collectively achieve Bose-Einstein condensation [2]. Their exceptionally light mass and the finite photonic lifetime make polaritons intrinsically driven-dissipative: continuous non-resonant pumping feeds an incoherent reservoir that relaxes into a phase-coherent condensate while radiative decay converts the condensate field into a measurable optical signal [3]. Such features raise the question of whether these condensates can be used to emulate associated with non-equilibrium stochastic growth of interfaces universality classes – a key concept in critical phenomena, whereby disparate systems can exhibit the same scaling behavior near transitions and at large scales [4]. Two-dimensional arrays of exciton-polariton condensates have been employed as XY-simulators [5]. The phase-locking in an array of several tens of condensates can be achieved on a picosecond time-scale, which makes the polariton system fastest of all known analogue simulators.

It turns out that the unique combination of tunable gain, losses, interactions, and intrinsic noise in polariton condensates hosts dynamical universality classes that

Figure 1. Schematic illustration of a one-dimensional lattice of interacting polariton condensates created by an array of non-resonant pump spots (orange). Each spot generates a localized condensate whose tunnel coupling J to its nearest neighbors has both coherent (Jr) and dissipative (Ji) components. These coupling strengths can be precisely tuned by adjusting the geometry and power of the pump spots.

have no analogue in systems at thermal equilibrium, where detailed balance precludes such non-linear scaling behaviour [6].

A particularly important nonequilibrium fixed point is associated with the KPZ universality class originally introduced to describe classical surface growth [7]. Inter-

¹Abrikosov Center for Theoretical Physics, MIPT, Dolgoprudnyi, Moscow Region 141701, Russia ²Russian Quantum Center, Skolkovo, Moscow, 121205, Russia ³Hybrid Photonics Laboratory, Skolkovo Institute of Science and Technology, 121205 Moscow, Russia ⁴Science Institute, University of Iceland, Dunhagi 3, IS-107 Reykjavik, Iceland

^∗ andrewkudlis@gmail.com

ested readers may consult Ref. [8] for a thorough review of kinetic roughening and the KPZ universality class. In one dimension the KPZ equation on h(x, t), representing the height (or, more generally, the phase-like interface displacement) field, reads

$$\partial_t h = \nu \partial_x^2 h + \frac{\lambda}{2} (\partial_x h)^2 + \eta(x, t),$$ (1)

where ν is an effective surface tension, λ is the non-linear growth coefficient, and η(x, t) is a zero-mean Gaussian white noise. This equation exhibits universal exponents β = 1/3 (temporal growth) and χ = 1/2 (spatial roughness). Their exact values in 1D were first derived in Ref. [7]. In the next section, we define them by the scaling laws. Mapping the nonlinear Schrödinger equation onto a noisy Burgers fluid allows using the condensates of one-dimensional weakly interacting bosons for emulation of the KPZ dynamics [9, 10]. A comprehensive account of these theoretical developments is also given in [6]. Subsequent exact results for generic Hamiltonian chains [11] and driven spin models [12] confirmed the ubiquity of this scaling. Early experimental confirmation beyond Eden-type growth came from measurements on turbulent liquid-crystal interfaces [13]. Direct experimental evidence of KPZ exponents in a single quasi-1D polariton channel was reported by Fontaine et al. [14]. Temporal coherence deviations from linear theory, matching KPZ scaling, are shown in Ref. [15]. KPZ scaling was also investigated in optical parametric oscillator regime of microcavity polaritons by Ferrier et al. [16]. More recently, Deligiannis et al. demonstrated that discretised two-dimensional polariton lattices also exhibit KPZ roughening in the presence of topological defects [17, 18]. Scaling properties of one-dimensional driven-dissipative condensates were also numerically studied in Ref. [19]. These results raise a fundamental and a practical question: under which conditions does a lattice of multiple interacting condensates support KPZ scaling in its phase field? Besides its conceptual interest, answering this question is crucial for analogue simulation schemes that rely on coherent networks of polariton condensates. In practice, the KPZ or EW character of the phase dynamics sets the attainable coherence length and computation time of a polariton-based analogue simulator; our phase diagrams provide a concrete recipe for choosing pump geometry and inter-site coupling so that the device operates either in a highly coherent (EW) regime or, when desired, in a noise- dominated (KPZ) regime for stochastic sampling tasks.

Moroney and Eastham [20, 21] considered a one-dimensional chain of exciton-polariton condensates with complex nearest-neighbour coupling J ≡ J^r + iJⁱ and showed that the imaginary component, which stems from reservoir-mediated dissipative tunnelling, can destabilise the synchronous phase via the proliferation of space-time vortices [20]. Their analytical bound suggests that KPZ behaviour should survive provided that |Jⁱ | is sufficiently small compared with Jr. However, a systematic exploration of the (Jr, Ji) phase diagram is missing in their study.

In the present work we derive the effective phase dynamics of an optically induced polariton lattice starting from a microscopic driven-dissipative Gross-Pitaevskii equation (GPE) that includes stochastic reservoir kinetics. In Fig. 1, we schematically present the model system. It represents a one-dimensional array of exciton-polariton condensates created by a non-resonant optical pumping. The condensates lose polaritons due to their radiative decay. They are replenished due to the stimulated scattering of excitons from incoherent reservoirs. We account to the coherent and dissipative couplings of the condensates [22, 23]. By adiabatically eliminating density fluctuations we obtain a stochastic nonlinear equation for the lattice phase which reduces to the KPZ form within a finite wedge of coupling parameters. We corroborate this analytical criterion by extensive numerical simulations of chains with up to 2000 condensates – matching present experimental capabilities – and analyse finite-size effects and the crossover to Edwards-Wilkinson scaling.

The remainder of the paper is organised as follows. Section II presents the driven-dissipative GPE with reservoir, its reduction to a discrete lattice model and the derivation of the KPZ coefficients. There, we also detail the numerical methods used in our simulations. Section III discusses the resulting phase diagrams and scaling exponents. Section IV summarises our findings and outlines perspectives of extension of the model for description of KPZ-simulators based on higher-dimensional lattices of polariton condensates.

II. MICROSCOPIC MODEL AND KPZ REDUCTION

In this section we formulate the driven-dissipative Gross-Pitaevskii description appropriate for a one-dimensional array of optically trapped excitonpolariton condensates, derive its discrete lattice representation, and outline the continuum limit leading to Kardar-Parisi-Zhang (KPZ) scaling.

A. Driven-dissipative Gross-Pitaevskii equation

We consider a high-quality planar semiconductor microcavity pumped non-resonantly by a mask-shaped laser profile that generates a set of N spatially separated exciton reservoirs. The coherent polariton order parameter Ψ(r, t) interacts with the incoherent reservoir density nR(r, t) and obeys the open-dissipative Gross-Pitaevskii equation (ODGPE) [24]

$$i\hbar \partial_t \Psi = \left[ -\frac{\hbar^2 \nabla^2}{2m} + \frac{i\hbar}{2} R n_R - \frac{i\hbar}{2} \gamma_C |\Psi|^2 + \hbar \alpha |\Psi|^2 \right] \Psi + \hbar \xi,$$ (2a)

$$\partial_t n_R = P(\mathbf{r}) - (\gamma_R + R|\Psi|^2)n_R.$$ (2b)

Here m is the lower-polariton effective mass, α is the polariton-polariton interaction constant, γ^C (γR) is the decay rate of the condensate (reservoir), R is the stimulated-scattering rate from the reservoir into the condensate, and P(r) denotes the spatially patterned nonresonant pump. Fluctuations are represented by complex Gaussian white noise ξ with correlator ⟨ξ ∗ (r, t)ξ(r ′ , t′ )⟩ = 2D⁰ δ(r − r ′ )δ(t − t ′ ).

When the reservoir relaxes much more slowly than the condensate, γ^R ≪ γ^C , the characteristic reservoir timescale τ^R = γ −1 ^R greatly exceeds the condensate timescale τ^C = γ −1 C . In this adiabatic limit n^R can be treated as quasi-static: nR(r, t) ≃ n (0) ^R (r). Dropping Eq. (2b) and replacing the gain term ⁱ^ℏ ² Rn^R in Eq. (3) by an effective pump P^C (r) = Rn(0) ^R , one obtains the reduced form

$$i\hbar \,\partial_t \Psi = \left[ -\frac{\hbar^2 \nabla^2}{2m} + \frac{i\hbar}{2} P_C - \frac{i\hbar}{2} \gamma_C |\Psi|^2 + \alpha |\Psi|^2 \right] \Psi + \xi, \quad (3)$$

This simplified equation is used in the subsequent derivation of the discrete lattice model.

B. Discrete lattice model

The patterned pump is arranged as an array of N identical Gaussian spots centred at positions xⁿ = na (n = 1, . . . , N) with lattice period a. Each spot supports a normalised single-particle orbital ϕn(r) whose spatial overlap is appreciable only with its nearest neighbours. Representing the condensate field as

$$\Psi(\mathbf{r},t) = \sum_{n=1}^{N} \psi_n(t) \,\phi_n(\mathbf{r}), \tag{4}$$

and projecting the reduced ODGPE (3) onto ϕⁿ yields a discrete driven-dissipative Gross-Pitaevskii equation for the complex amplitudes ψn(t).

After retaining nearest-neighbour hopping and adopting the notation introduced earlier, the evolution equation that is actually integrated reads [20]

$$i\hbar\dot{\psi}{n} = \left[E{0} - \frac{i\hbar}{2}\gamma{0}|\psi{n}|^{2} + \frac{i\hbar}{2}P + \alpha|\psi{n}|^{2}\right]\psi{n} + \sum{k=\langle n\rangle}J\psi{k} + \hbar\xi_{n}, \quad (5)$$

where E⁰ is the on-site kinetic energy, γ⁰ the radiative decay rate, P the effective linear gain supplied by the pump, α the on-site interaction strength, and J = Jr+iJⁱ the complex nearest-neighbour coupling introduced in Sec. II C. The stochastic term ξn(t) is the lattice counterpart of the white noise ξ(r, t) defined in Sec. II A. Equation (5) is solved numerically for chains containing up to N = 2000 condensates.

C. Connection to KPZ universality

All numerical data reported below are obtained by integrating the complex-field lattice equation (5). To interpret those results it is instructive to figure out how the same microscopic model reduces, at long wavelengths, to the one-dimensional Kardar-Parisi-Zhang (KPZ) equation that governs the universal fluctuations of the condensate phase.

We decompose the lattice field into amplitude and phase

$$\psi_i(t) = \sqrt{n_i(t)} e^{i\theta_i(t)}, \tag{6}$$

where nⁱ and θⁱ are the on-site density and phase. Density fluctuations relax on the time-scale of τⁿ ∼ (γ0) −1 , which is short compared to the time-scale of phase diffusion. Eliminating δnⁱ = nⁱ −n adiabatically and expanding finite differences to the second order in the lattice spacing a one arrives, in the continuum limit x = ia, at

$$\partial_t \theta(x,t) = \nu \,\partial_x^2 \theta + \frac{\lambda}{2} (\partial_x \theta)^2 + \zeta(x,t),$$ (7)

the standard one-dimensional KPZ equation. The effective diffusion and non-linear coefficients read

$$\nu = \frac{2a^2(J{\rm i} - \eta J{\rm r})}{\hbar}, \qquad \lambda = -\frac{a^2(J{\rm r} + \eta J{\rm i})}{\hbar}, \qquad (8)$$

with η = 2α/(ℏγ0) and J = J^r + iJⁱ . The term ζ(x, t) is Gaussian white noise, inherited from ξi(t).

Equation (7) belongs to the strong-coupling KPZ universality class provided the diffusion constant is positive,

$$J{\rm i} > \eta J{\rm r}. \tag{9}$$

When, in addition, the non-linear coefficient λ vanishes (J^r +ηJⁱ = 0) the dynamics crosses over to the Edwards-Wilkinson universality class with critical exponents β = 1/4 and χ = 1/2. Although this coarse-grained reduction is not used in the simulations, condition (9) provides a practical guideline for selecting hopping parameters that display KPZ scaling in the full complex-field model.

D. Numerical integration protocol

To capture nonlinear effects beyond a phase-only description we integrate the full complex-field equation (5). In our simulations, the chain length is varied up to N = 2000, enabling a systematic finite-size study.

The stochastic differential equations are treated with a semi-implicit Euler-Maruyama algorithm. Space-time first-order coherence is evaluated on the lattice

$$g^{(1)}(\Delta x, \Delta t) = \frac{\langle \psi^*(x_0, t_0) \psi(x', t') \rangle}{\langle |\psi(x_0, t_0)| |\psi(x', t')| \rangle}, \tag{10}$$

where x ′ = x⁰ + ∆x, t ′ = t⁰ + ∆t, and ⟨· · · ⟩ denotes an ensemble average. The the coordinate-time pair (x0, t0) is chosen deeply inside the chain after transients have decayed.

To improve the ensemble statistics, we calculate the ensemble average over M different systems of N condensates by considering chain of condensates with length N · M, which consists of M independent sub-chains of the length N each. Time evolution and average values of this ensemble were calculated on a GPU using Compute Unified Device Architecture (CUDA) library.

The critical exponents are then extracted from the following scaling relations

$$C(\Delta x = 0, \Delta t) \equiv -2\ln|g^{(1)}(0, \Delta t)| \propto \Delta t^{2\beta}, \qquad (11)$$

$$C(\Delta x, \Delta t = 0) \equiv -2 \ln |g^{(1)}(\Delta x, 0)| \propto \Delta x^{2\chi},$$ (12)

with the exponents β and χ obtained by least-squares fits over two decades in ∆t or ∆x on log-log axes.

Throughout we set ℏ = 1 and measure energy in units of the polariton radiative decay rate γ0; physical time is therefore rescaled as τ = tγ0. A dimensionless time step of ∆τ = 0.1 was adopted, and convergence was explicitly verified by repeating the simulations with smaller steps. The linear pump is fixed at P/γ⁰ = 1, and the interaction-to-damping ratio is kept weak, η = 2α/(ℏγ0) = 10⁻⁴ , as appropriate for largearea microcavities. Also note, that the on-site kinetic energy E⁰ does not affect first-order correlator g (1) and only adds a constant drift in a phase temporal evolution (θi(t) ∼ −E0t/ℏ), and is therefore set to zero without loss of generality. Within this parameter set we scan the complex coupling amplitude J in both real and complex directions. This variation can be achieved by changing the pump geometry. Experimental engineering of spatial coherence and couplings in polariton lattices appears in Ref. [25].

The noise strength D can be tuned to change the coherence time of the system but it should be low enough for condensate fluctuation being small (∆nⁱ ≪ ni). This condition is required by the adiabatic elimination that leads to the KPZ equation (7). We choose D so that the root-mean-square density fluctuation equals 5 % of the steady-state density obtained in the noiseless limit

$$\sqrt{D} = 0.05 \frac{P + 4J_{\rm i}}{\sqrt{2\gamma_0 \Delta \tau}}.$$ (13)

An initial condition for ψ does not affect steady state solution at t → ∞ but it defines how fast the system's dynamics enters the specific universality class regime. We are using equal-phase initial conditions ψi(0) = ψ^j (0) = const. With such initial conditions the system is spatially correlated at t = 0: g (1)(∆x, 0) ^t0=0 = 1. So KPZ universality class characterized by spatial correlations of Eq.(12) can be established only after some time of phase fluctuations. All statistical sampling begins after a waiting time of τ⁰ = 5 × 10⁴ , which is well beyond this transient.

III. RESULTS

In this section, we explore how the critical exponents extracted from Eq. (5) depend on the complex hopping amplitude J = J^r +iJⁱ , on the interaction-damping ratio η = 2α/(ℏγ0), and on the chain length N. Unless stated otherwise the interaction is fixed at the experimentally relevant value η = 10−⁴ and the noise strength is chosen as described in Sec. II D.

First, as argued at the end of Sec. II D, the KPZ scaling regime is reached only after an initial transient. Before running the large-scale simulations we therefore determine a waiting time t⁰ beyond which correlations can

Figure 2. Relaxation of the spatial coherence toward the KPZ regime. Panel (a): equal-time correlator g (1)(∆x, 0) computed for different waiting times τ⁰ = t0γ0: 0, 2×10⁴ , 4×10⁴ , 6 × 10⁴ , 8 × 10⁴ , and 10⁵ . Panel (b): roughness exponents χ extracted from the large-∆x tails. For τ⁰ = 5×10⁴ the correlator has completely lost memory of the uniform initial state and the exponent has approximately converged to the KPZ value χ ≃ 1/2 (horizontal black dashed line), justifying the choice of this waiting time for all subsequent analysis. The selection of τ⁰ is crucial to ensure that transient effects from the initial coherent phase are fully relaxed, so that subsequent measurements reflect the true KPZ scaling regime.

Figure 3. (a) Colour map of the temporal growth exponent β extracted from −2 ln |g (1)(0, ∆t)| ∝ ∆t 2β in the plane of real and imaginary parts of the hopping amplitude (Jr, Ji). The data correspond to interaction η = 10⁻⁴ , chain length N = 500, and ensemble size M = 10 000. (b) The same data represented in the variables λ = −(Jr+ηJi) and ν = 2(Ji−ηJr) that enter the coarse-grained KPZ equation (recall that ℏ = 1 and a = 1). Because the working value η= 10⁻⁴ ≪ 1, the mapping produces, to excellent accuracy, a mirror reflection of panel (a) with respect to the vertical axis; we include it to make the KPZ phase boundary immediately visible in the conventional (λ, ν) language. (c) Log-log plots of −2 ln |g (1)(0, ∆t)| for the two marked points A (light blue pint on panel (a) and (b)) (J/γ⁰ = −0.2+ 0.03i), where KPZ scaling is well pronounced, and B (deep blue point on both upper panels) (J/γ⁰ = 0.04i), for which EW scaling is expected. Linear fits (dashed) yield β^A ≃ 0.325 (1/3) for A, characteristic of the KPZ class, and β^B ≃ 0.248 (1/4) for B, consistent with Edwards-Wilkinson (EW) scaling. Panel (c) was obtained for a larger ensemble (M = 100 000) to suppress statistical noise.

be analysed safely. Fig. 2 shows the equal-time spatial correlator g (1)(∆x, 0) and the roughness exponent χ extracted from it for several choices of t0. The data demonstrate that the value adopted in the following, τ⁰ = t0γ⁰ = 5 × 10⁴ , is already sufficient for the influ-

Figure 4. Spatial roughness exponent χ. (a) Colour map of χ in the (Jr, Ji) plane for interaction η = 10⁻⁴ , chain length N = 500, and M = 10 000 noise realizations. (b) Same data expressed in the KPZ variables λ and ν; because η ≪ 1 the image is essentially a mirror reflection of panel (a). (c) Log–log plot of −2 ln |g (1)(∆x, 0)| for points A (J/γ⁰ = −0.2 + 0.03i) and B (J/γ⁰ = 0.04i), where a linear fit (dashed lines) gives χ^A ≃ 0.528 ≈ 1/2 and χ^B ≃ 0.473 ≈ 1/2, in agreement with KPZ and EW scaling.

ence of the uniform initial condition to vanish. Therefore, this time moment is used as the starting point for all statistical analyses.

Now let us emphasize the main result of our work. Fig. 3(a) shows the colour map of the temporal growth exponent β in the complex J plane. For the chosen weak interaction the KPZ criterion Jⁱ > ηJ^r (see sec. II C) nearly corresponds to the upper half-plane Jⁱ > 0, except for the vertical axis J^r = 0, where the non-linear coefficient λ vanishes and Edwards-Wilkinson (EW) scaling should be observed. This analytical expectation is confirmed by the simulations: a narrow red stripe with β ≃1/4 (EW class) appears along J^r = 0 axis, while the bulk of the upper half-plane is violet with β ≃1/3 (KPZ class). The green area represents a crossover region: the

Figure 5. Dependencies of the critical exponents on the size of the system at J/γ⁰ = −0.2 + 0.03i (point A in Fig. 3 and Fig. 4). Top panel: temporal growth exponent β as a function of the chain length N; the horizontal dashed line marks the KPZ value for exponent βKPZ = 1/3, while the solid orange line shows the numerical plateau β^N = 0.329±0.003 obtained by averaging the last twenty points (N ≥1000). Bottom panel: spatial roughness exponent χ versus N. The dashed line indicates the KPZ value χKPZ = 1/2 and the solid orange line shows the numerical estimate χ^N = 0.504±0.004. Each point represents an average over M = 10⁴ independent noise realizations.

fitted exponents there differ from both KPZ and EW because the required crossover scales exceed our simulation window. In the infinite-size, long-time limit these points are expected to flow to the EW fixed point. Close to the real axis the fit of −2 ln |g (1)(0, ∆t)| ∝ ∆t ²^β becomes unreliable (light blue region) because the diffusive coefficient ν is very small, and the statistical spread of the phase grows, requiring much more noise realizations for convergence. Fig. 3(b) replots exactly the same data in (λ, ν) variables of the KPZ equation using the relations λ = −(J^r + ηJi) and ν = 2(Jⁱ − ηJr), where we allow a = 1. With our experimentally realistic value η = 10⁻⁴ the transformation reduces essentially to a mirror reflection with respect to the vertical axis, so panel (b) is the left–right image of panel (a); nevertheless it makes evident the location of the KPZ wedge in the (λ, ν) plane. Fig. 3(c) illustrates the two universality classes explicitly. For the point J/γ⁰ = −0.2 + 0.03i the log − log plot of −2 ln |g (1)(0, ∆t)| yields a slope 2β = 2/3, characteristic of KPZ scaling, whereas for the purely imaginary coupling J/γ⁰ = 0.04i the slope is 2β = 1/2, as predicted by EW theory.

Fig. 4 shows the counterpart of Fig. 3(a) for the spa-

tial roughness exponent χ, extracted from the scaling law −2 ln |g (1)(∆x, 0)| ∝ ∆x ²χ. Panel (a) shows the colour map in the original variables (Jr, Ji). Over nearly the entire upper half-plane (Jⁱ > 0) the simulations predict χ ≃ 0.5, in the excellent agreement with the onedimensional KPZ value χKPZ = 1/2. Close to the real axis, however, a narrow red band appears where χ cannot be determined reliably; this is the same region where Fig. 3(a) showed large uncertainties in the exponential factor β. In the limit Jⁱ → 0 the effective diffusion constant ν vanishes, phase gradients become extremely rough, and the number of noise realizations required for converged statistics grows faster than our ensemble size. The light blue EW region at J^r = 0 in Fig. 4 likewise originates from the finite system size. We note that for the spatial correlator the roughness factor χ takes the same universal value in both the KPZ and EW classes. When the real part of the coupling is set to zero the numerical convergence deteriorates and finite–size effects of the simulated chain become much more pronounced. For completeness we re-plot exactly the same data in panel (b) using the KPZ parameters λ = −(J^r +ηJi) and ν = 2(Jⁱ − ηJr). Because the working value η = 10⁻⁴ is extremely small, it amounts practically to a mirror reflection of panel (a) about the vertical axis. This alternative representation makes the location of the KPZ wedge in the conventional (λ, ν) plane quite clear. Panel (c) illustrates the scaling at the representative point A (marked in panels a,b). The log–log plot of −2 ln |g (1)(∆x, 0)| versus ∆x/a shows a slope 2χ ≃ 0.966, i.e. χ ≃ 0.483, that is close to the KPZ value given the numerical accuracy imposed by the size of the finite system.

Fig. 5 (top panel) tracks the growth exponent β as a function of chain length for the representative coupling J/γ⁰ = −0.2 + 0.03i. Starting from β ≈ 0.31 at N = 100, the value rises steadily and saturates at the KPZ benchmark βKPZ = 1/3 once N ≳ 1000. The bottom panel of Fig. 5 shows the companion finite-size flow of χ, which approaches 1/2 on the same length scale. These trends confirm that the EW-like window shrinks as N ⁻¹ and collapses onto the line J^r = 0 in the thermodynamic limit, as it is expected from the analytical criterion Jⁱ > ηJr.

We note, however, that a residual spread of individual data points is still visible in both panels. This spread lies within the numerical accuracy of the simulation and reflects (i) the finite number of noise realizations, (ii) the stochastic nature of the least–squares fits, and (iii) the intrinsic difficulty of extracting scaling exponents from a highly nonlinear driven–dissipative system.

To provide definite numerical benchmarks for the most favorable parameter region highlighted in Figs. 3 (c) and 4 (c), we quote final estimates obtained as the average over the last twenty data points in Fig. 5, i.e. over chains with N ≥ 1000 where additional condensates do not improve the statistics further. The resulting values

$$\beta{N,\text{KPZ}} = 0.329(3), \qquad \chi{N,\text{KPZ}} = 0.504(4)$$ (14)

are drawn as solid orange guide lines in the figure and are taken as the the final numerical KPZ estimates of our model.

In summary, both critical exponents approach their KPZ benchmarks once the chain hosts only a few hundred condensates – an array size that can be realized with current spatial–light–modulator techniques and state-of-theart microcavities. The convergence is monotonic: as N increases the growth exponent β rises from its finite-size value and locks onto 1/3, while the roughness exponent χ descends toward 1/2 on the same length scale, indicating that phase gradients have reached their asymptotic strong-coupling statistics. The modest residual offsets observed very close to Jⁱ = 0 should not be interpreted as a failure of KPZ scaling. They stem from the vanishing diffusion constant ν ∝ Jⁱ , which slows down the equilibration of long-wavelength modes and amplifies statistical noise. Longer simulation times or a larger ensemble of noise realizations would reduce this bias further, but the overall trends remain unchanged.

IV. CONCLUSION

We have combined analytical and numerical approaches to chart the conditions under which onedimensional arrays of driven–dissipative excitonpolariton condensates emulate the Kardar–Parisi–Zhang universality class equations. Starting from a stochastic Gross-Pitaevskii description that retains both gain saturation and reservoir noise, we show that adiabatic elimination of density fluctuations yields an effective KPZ equation for the lattice phase whenever the inequality Jⁱ > ηJ^r is satisfied. Extensive GPU simulations – covering chain lengths up to two thousand sites and ensembles of 10⁴ noise realisations – corroborate this criterion: the growth and roughness exponents approach the KPZ benchmarks, β = 1/3 and χ = 1/2, once the system contains a few hundred condensates. For the coupling J/γ⁰ = −0.2 + 0.03i we obtain final finite-size estimates of β^N = 0.329(3) and χ^N = 0.504(4), providing the specific numerical targets for experiments. Our findings reconcile the KPZ scaling observed in single polariton wires [14] with the space–time-vortex scenario proposed for coupled condensates [20], and delineate the parameter window where coherent polariton networks can serve as analogue simulators of non-equilibrium surface growth. Further studies are required to the present analysis to two-dimensional lattices, explore the impact of disorder and detuning, and investigate whether engineered pump geometries can stabilise higher-dimensional KPZ fixed points in driven-dissipative quantum fluids.

Acknowledgments.– The research is supported by the Ministry of Science and Higher Education of the Russian Federation (Goszadaniye) Project No. FSMG-2023-0011. The work of A.K. is supported by the Icelandic Research Fund (Rannsóknasjóður, Grant No. 2410550).

[1] I. Carusotto and C. Ciuti, Quantum fluids of light, Rev. Mod. Phys. 85, 299 (2013).
[2] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose–Einstein condensation of exciton polaritons, Nature 443, 409 (2006).
[3] P. B. Littlewood and A. Edelman, Universal Themes of Bose–Einstein Condensation, edited by N. P. Proukakis, D. W. Snoke, and P. B. Littlewood (Cambridge University Press, Cambridge, 2017).
[4] H. Stanley, Phase Transitions and Critical Phenomena (Oxford University Press, 1971).
[5] K. P. Kalinin and N. G. Berloff, Simulating Ising and n-state planar Potts models and external fields with nonequilibrium condensates, Phys. Rev. Lett. 121, 235302 (2018).
[6] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Rep. Prog. Phys. 79, 096001 (2016).
[7] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56, 889 (1986).
[8] T. Halpin-Healy and Y.-C. Zhang, Kinetic roughening phenomena, stochastic growth, directed polymers and all that. aspects of multidisciplinary statistical mechanics, Phys. Rep. 254, 215 (1995).
[9] M. Kulkarni and A. Lamacraft, Finite-temperature dynamical structure factor of the one-dimensional Bose
gas: From the Gross–Pitaevskii equation to the Kardar– Parisi–Zhang universality class of dynamical critical phenomena, Phys. Rev. A 88, 021603(R) (2013).
[10] J. De Nardis, S. Gopalakrishnan, and R. Vasseur, Nonlinear fluctuating hydrodynamics for Kardar-Parisi-Zhang scaling in isotropic spin chains, Phys. Rev. Lett. 131, 197102 (2023).
[11] H. van Beijeren, Exact results for anomalous transport in one-dimensional hamiltonian systems, Phys. Rev. Lett. 108, 180601 (2012).
[12] B. Derrida and C. Appert, Universal large-deviation function of the Kardar–Parisi–Zhang equation in one dimension, J. Stat. Phys. 94, 1 (1999).
[13] K. A. Takeuchi, M. Sano, T. Sasamoto, and H. Spohn, Growing interfaces uncover universal fluctuations behind scale invariance, Sci. Rep. 1 (2011).
[14] Q. Fontaine, D. Squizzato, F. Baboux, I. Amelio, A. Lemaître, M. Morassi, I. Sagnes, L. Le Gratiet, A. Harouri, M. Wouters, I. Carusotto, A. Amo, M. Richard, A. Minguzzi, L. Canet, S. Ravets, and J. Bloch, Kardar–Parisi–Zhang universality in a onedimensional polariton condensate, Nature 608, 687 (2022).
[15] K. Ji, V. N. Gladilin, and M. Wouters, Temporal coherence of one-dimensional nonequilibrium quantum fluids, Phys. Rev. B 91, 045301 (2015).
[16] A. Ferrier, A. Zamora, G. Dagvadorj, and M. H. Szymańska, Searching for the Kardar-Parisi-Zhang phase in
microcavity polaritons, Phys. Rev. B 105, 205301 (2022). [17] K. Deligiannis, Q. Fontaine, D. Squizzato, M. Richard, S. Ravets, J. Bloch, A. Minguzzi, and L. Canet, Kardar-Parisi-Zhang universality in discrete two-dimensional driven-dissipative exciton polariton condensates, Phys. Rev. Res. 4, 043207 (2022).
[18] K. Deligiannis, D. Squizzato, A. Minguzzi, and L. Canet, Accessing Kardar-Parisi-Zhang universality sub-classes with exciton polaritons(a), Europhys. Lett. 132, 67004 (2021).
[19] L. He, L. M. Sieberer, E. Altman, and S. Diehl, Scaling properties of one-dimensional driven-dissipative condensates, Phys. Rev. B 92, 155307 (2015).
[20] J. P. Moroney and P. R. Eastham, Synchronization and spacetime vortices in one-dimensional driven-dissipative condensates and coupled oscillator models, Phys. Rev. B 108, 195302 (2023).
[21] J. P. Moroney and P. R. Eastham, Synchronization in disordered oscillator lattices: Nonequilibrium phase transition for driven-dissipative bosons, Phys. Rev. Res. 3, 043092 (2021).
[22] I. L. Aleiner, B. L. Altshuler, and Y. G. Rubo, Radiative coupling and weak lasing of exciton-polariton condensates, Phys. Rev. B 85, 121301 (2012).
[23] H. Ohadi, R. L. Gregory, T. Freegarde, Y. G. Rubo, A. V. Kavokin, N. G. Berloff, and P. G. Lagoudakis, Nontrivial phase coupling in polariton multiplets, Phys. Rev. X 6, 031032 (2016).
[24] M. Wouters and I. Carusotto, Excitations in a nonequilibrium Bose-Einstein condensate of exciton polaritons, Phys. Rev. Lett. 99, 140402 (2007).
[25] J. D. Töpfer, I. Chatzopoulos, H. Sigurdsson, T. Cookson, Y. G. Rubo, and P. G. Lagoudakis, Engineering spatial coherence in lattices of polariton condensates, Optica 8, 106 (2021).

novokreschenovkardarparisizhanguniversality2025.md 39 KB History Raw