Synchronization in disordered oscillatory media: a nonequilibrium phase transition for driven-dissipative bosons

We show that a lattice of phase oscillators with random natural frequencies, described by a generalization of the nearest-neighbor Kuramoto model with an additional cosine coupling term, undergoes a phase transition from a desynchronized to a synchronized state. This model may be derived from the complex Ginzburg-Landau equations describing a disordered lattice of driven-dissipative Bose-Einstein condensates of exciton polaritons. We derive phase diagrams that classify the desynchronized and synchronized states that exist in both one and two dimensions. This is achieved by outlining the connection of the oscillator model to the quantum description of localization of a particle in a random potential through a mapping to a modified Kardar-Parisi-Zhang equation. Our results indicate that long-range order in polariton condensates, and other systems of coupled oscillators, is not destroyed by randomness in their natural frequencies.

Introduction.-Synchronization of coupled oscillators is a phenomenon that appears regularly throughout nature. Many seemingly unrelated systems which exhibit repetitive behaviors, such as clocks, pacemaker cells in the heart, or a swarm of pulsing fireflies, are seen to undergo transitions from initial randomness to an ordered state [1,2]. A celebrated model of such synchronization phenomena was introduced by Kuramoto [3,4]. With allto-all couplings the Kuramoto model undergoes a phase transition, from a desynchronized state at weak coupling to a synchronized state at strong coupling. However, such ordering is absent when the oscillators are on a one or two dimensional lattice, with coupling between neighboring sites [5][6][7].
Synchronization plays an important role in the physics of Bose-Einstein condensation in driven-dissipative Bose gases. Such condensation has been realized for exciton polaritons [8], which are bosonic quasiparticles formed from the strong coupling of excitons and photons in semiconductor microcavities. The condensates are described by a single macroscopic wavefunction [9], giving rise to phenomena such as superfluidity, Josephson oscillations [10], and quantized vortices, and enabling applications such as analog simulation [11,12]. In these systems the decay of the polaritons is offset by gain from the pump, and condensation occurs in a nonequilibrium steady state. Above threshold the nonlinear gain fixes the density of the condensate, forming an autonomous phase oscillator whose frequency corresponds to the condensate's energy. Spatially separated condensates can form in the random potentials arising from disorder in the samples [13], and in engineered potentials such as lattices [11,14]. If the condensates are well separated they oscillate independently, but their mutual coupling can lead them to synchronize to a common frequency and phase [15]. Such effects have been observed for polaritons in double-well potentials [10,15], in the intrinsic random potential of the samples [13,16], and in a weakly disordered two dimensional lattice [14].
In this Letter we analyze the coherence properties of driven-dissipative bosons in a disordered lattice potential, focusing on the two dimensional case realized in experiment [14]. We seek to understand the phase diagram and phase transitions, as generalizations of the equilibrium case, in which there is a superfluid-insulator transition [17]. The nonequilibrium problem is described by a lattice of coupled oscillators, with disordered frequencies, suggesting that this generalization may be a frequencyordering transition, or a frequency-and phase-ordering transition. Ordered states are, however, not expected for Kuramoto oscillators with random frequencies in two dimensions [5,7]. This is consistent with a perturbative treatment of the effect of disorder on the polariton condensate, which predicts the absence of long-range phase order [18].
Here, we investigate synchronization in a lattice of driven-dissipative condensates with random natural frequencies. Such a lattice is described by a generalization of the nearest neighbor Kuramoto model that includes an additional coupling term that is even in the relative phases. Such models have recently been considered [19,20] in the case where the natural frequencies are identical and frequency synchronization is expected, but the phases may become disordered due to the spaceand time-dependent noise associated with gain and loss [21][22][23][24]. We consider the opposite limit, of random natural frequencies but negligible time-dependent noise. We find that the non-Kuramoto coupling has a significant effect: it leads to a true synchronization transition to a state with long-range frequency and phase order for dimensions d < 4. Our argument involves relating the oscillator model to a Kardar-Parisi-Zhang (KPZ) equation [25] with time-independent noise [26][27][28][29][30], and then to an imaginary-time Schrödinger equation with a random potential. This allows us to derive the phase boundary for synchronization and characterize the frequency and phase profiles. These analytical predictions agree with numerical simulations in one and two dimensions. Our results show that there is a nonequilibrium generalization of the superfluid-insulator transition in a lattice of driven-dissipative condensates, and that there is a region where there is long-range coherence that is robust against disorder. Our conclusions apply more generally to coupled oscillator systems, implying there are many other settings [1] in which this phase transition could be realized.
Dynamics of nonequilibrium condensates.-Under the requisite conditions [31], a polariton condensate may be described by a driven-dissipative Gross-Pitaevskii equation [32], where U 0 is the polariton-polariton interaction strength, while g 0 (r) and Γ 0 account for the linear gain -resulting from the difference of pumping and decay -and the gain saturation. V 0 (r) is a confining potential, which can arise from the repulsive interaction with the exciton reservoir, etching, deposition, and the intrinsic disorder in the sample [14,33]. We seth = 1 throughout. We consider a lattice of N condensates, as realized experimentally [14]. In that experiment a spatially patterned pump beam populates the exciton reservoir, forming a confining lattice potential. This potential is superimposed on the random potential due to disorder in the sample, leading to small variations in the energy from site to site. Condensation occurs in the lowest energy state of each well of the lattice potential. We model this by expanding the macroscopic wavefunction over the basis set of wavefunctions localized in individual wells [34,35]. Assuming that the overlap between neighboring wavefunctions is small results in N coupled equations for the amplitudes of each site, ψ k , Here, k is the energy of the condensate on a given lattice site, U , g and Γ are the interaction strength and gain coefficients, and the sum is over nearest neighbors. J kl ≈ J is the matrix element describing tunneling between neighboring sites. Since the pump pattern, which produces the largest contribution to the potential and controls the gain on each site, is periodic, these parameters will vary little from site to site, and we treat them as constant. The site energies however vary randomly from site to site, with their standard deviation σ providing a measure of the strength of the disorder.
To examine synchronization in this system, we reparameterize our equations in terms of density and phase , ψ k = √ n k exp (−iθ k ). Well above threshold, the density variation between condensates is small, δn kl = n k − n l n k . Furthermore, the fast relaxation of the densities allows for their adiabatic elimination. This gives [35] (3) Here, we have introduced the dimensionless parameter α ≡ Γ/U . The overall blueshift g/α can be removed by a redefinition of the zero of frequency. We can choose 1/σ as our unit of time, so that the solutions to Eq. (3) are controlled by two dimensionless parameters, α and J/σ. Eq. (3) is an equation for the phase dynamics in a system of coupled self sustained oscillators. If we neglect the cosine term in the sum, it is the nearest-neighbor Kuramoto model [5]. In general, we expect that the existence of a synchronized solution -whereby the oscillators rotate at a common frequency,θ 1 =θ 2 = · · · =θ N = Ω -will be dependent on the magnitudes of the tunneling and the spread of on-site energies. If there is no tunneling (J = 0), each oscillator will rotate at its blueshifted natural frequency. When J is large relative to the spread of the natural frequencies however, the frequency of each oscillator is strongly affected by the phase difference between it and its neighbors, and this mechanism can bring about synchronization. Despite this, it has been shown that one and two dimensional lattices of Kuramoto oscillators, with random on-site energies and nearest-neighbor couplings, do not exhibit synchronization [5][6][7] in the limit N → ∞. In one dimension the probability of synchronization can be calculated as a function of N , and a coupling strength of order √ N is required for a synchronized solution [6].
The presence of the cosine term in the coupling function of our model has a significant effect however. Unlike the Kuramoto model, our coupling is non-odd in its ar- guments, and it has been suggested that this may bring about synchronization more readily [36,37]. Numerical simulations of Eq. (3) illustrate that this is indeed the case. Fig. 1 shows the probability of synchronization for chains of condensates of various lengths with normallydistributed on-site energies. This probability is determined by calculating the average frequency of each oscillator from its phases at two times, one a short time after the initial transient behavior has decayed, and the other a long time later. If each bond between neighboring sites has a frequency difference less than the smallest numerically resolvable frequency the configuration is considered synchronized. P sync (J) is then estimated from the fraction of synchronized results arising over many disorder realizations. We see that the probability of synchronization varies from almost zero to almost one over a range of J. The width of this range is non-zero, because different realizations synchronize at slightly different tunneling strengths. More importantly, the position of the center of this range, which we use to define a typical critical tunnneling strength J c such that P sync (J c ) = 0.5, is almost independent of the system size. This strongly suggests that there is a synchronization transition in the limit of thermodynamically large systems. This would be a non-equilibrium phase transition, at which the steadystate of a system in the thermodynamic limit changes character, giving rise to singularities in its properties. Such behavior is markedly different from that of the Kuramoto model, where J c scales with √ N , so that synchronization occurs only in small systems.
Continuum theory-To understand the synchronization transition, and the nature of the synchronized states, we consider the continuum form of Eq. (3) in the limit where the phase differences between neighboring sites are small. This will be the case in a synchronized state for weak disorder. Expanding the trigonometric functions to second order, and taking the continuum limit, we have where a is the lattice constant, which will be set to one in the following. A uniform energy shift g/α − 2Jd has been absorbed in the definition of θ(x, t), d being the dimensionality of the system -1 or 2 for the cases of interest to us. Eq. (4) is similar to the KPZ equation, but has a time independent noise term [25,29]. The connection between the complex Ginzburg-Landau equation and the conventional KPZ equation, with spatiotemporal noise, has been made in previous studies of polariton condensation [21][22][23][24], and indeed lattice models similar to Eq. (3) have also been studied [19,20,38]. Those works consider noise associated with gain and loss, rather than that due to a random potential. The form of the noise in Eq. (4) results from our consideration of purely spatial disorder, and leads to different universal behavior. The phase, θ(x, t), behaves like the height of an interface, with a growth rate (x) which is random in space but not in time.
A Cole-Hopf transformation enables us to write the KPZ equation as This is the imaginary time Schrödinger equation for a particle of mass α/2J in a random potential V (x) = −α (x). In this form it also describes the evolution of a population with diffusion and random autocatalytic amplification [26,27], and the partition function for a directed polymer in a random potential [28,30]. The general solution of Eq. (6) can be expressed in terms of the eigenstates ofĤ, with energies E n , which approaches the ground state wavefunction of the random potential V (x) for t → ∞. For long but finite times, Z(x, t) will have contributions from a small number of low-energy states. In dimensions d < 4 these lowlying states will be localized, ψ n ∼ e −|x−xn|/ζ , at some dilute positions x n , with the localization length ζ given by balancing the kinetic and potential terms in Eq. (6): [27]. Here D = J/α and u L = ασ, so we have .
In principle the ground state localization length should be obtained from the depth of the deepest of the ∼ (L/ζ) d potential wells of size ∼ ζ in the sample, rather than the depth of a typical well. These differ by a factor d ln L/ζ which, while it formally diverges as N = L d → ∞, is of order one in realistic systems [27]. In Figs. 2 and 3 we show some phase and frequency profiles obtained by solving the coupled oscillator model, Eq. (3). These solutions agree qualitatively with those of the continuum model, given by Eqs. (5) and (7). At short times, multiple peaks are evident in the phase profile, each of which corresponds to a localized, low-energy state of the Hamiltonian. The corresponding energies appear as plateaus in the frequency profile. The low-lying states have negative energies, and so they grow in magnitude, corresponding to a steadily increasing phase in the region controlled by each localization center. The ground state, with the lowest energy, grows fastest, and eventually dominates the solution, giving a solution with a single peak in the phase profile, and a single frequency (corresponding to the ground state energy ofĤ).
Phase diagram-The solutions to the continuum model always synchronize, given sufficient time, but this is not the case for the original model. To investigate the conditions for synchronization we consider the phase gradients. Since the solutions to Eq. (6) are formed from exponentially localized states with localization length ζ we have, from Eqs. (5) and (8), .
The continuum model does not impose a limit on these gradients, however, it only captures the behavior of the original model when |∇θ| < ∼ 1. Otherwise the periodicity and limited range of the trigonometric functions are important, and we cannot expand them in a power series. This argument gives as the phase boundary for synchronization. This condition is relevant for d < 4, as localized solutions to Eq. (6) are not guaranteed otherwise.
The dependence of J c on α given by Eq. (10) is confirmed by simulations for lattices in one and two dimensions, as shown in Fig. 4. We simulate Eq. (3) to obtain P sync (J), fit the resulting functions to a Gumbel extreme-value distribution, and use the fitted parameters to find J c . This extreme-value distribution is expected from the analysis above, since the system is synchronized at a given J if the ground-state localization length exceeds a certain value, or equivalently if the ground-state energy exceeds a certain value. It fits our results well, as can be seen in Fig. 1. The points in Fig. 4 show the values of J c . The solid curves show the predicted squareroot (1D) and linear (2D) dependencies, which are in good agreement with the data. On the 1D phase diagram, we also show the prediction of the conventional Kuramoto model for a chain of N = 800 sites, for which J c ∼ σα √ N . While the critical coupling in that case, without the cosine term, diverges as N → ∞, for a finite system it can nonetheless lie below Eq. (10) at small α. This explains the crossover seen at small α in Fig.  4a. This is the standard finite-size behavior if the cosine term is relevant in the renormalization group sense: such terms, even if they are small at the lattice scale, grow with distance, and control the physics in the limit N → ∞. In a finite system, however, the growth can be cut off by the system size.
Although our numerical results agree with Eq. (10) for α < ∼ 1, they disagree in the opposite regime, where we find large sample-to-sample fluctuations and many states which are desynchronized even at large J. This may be related to dynamical instabilities in that regime [19,20]. In any case, Eq. (8) holds only in the weak disorder regime where the localization length ζ is larger than the lattice spacing a = 1, and the localization length at the transition is, from Eqs. (8) and (10), ζ c ∼ 1/α. More generally, ζ c > 1 is needed so that space can be treated as continuous, as in Eq. (4), at the phase boundary. Importantly, the existence of a synchronization transition in the continuum regime implies the phenomenon is universal, and independent of the details of the lattice or disorder. We note that a synchronization transition for one dimensional chains of oscillators with non-odd nearestneighbor coupling was previously identified byÖstborn [37]. That analysis, however, predicts a critical J c which differs from Eq. (10), and does not agree well with our numerical results.

Conclusions-
We have shown that a model for coupled phase oscillators, which describes a disordered lattice of polariton condensates, has a synchronization transition in the thermodynamic limit. At this transition tunneling between condensates overcomes the localizing effects of the random potential, leading to a state with long-range coherence. This is reminiscent of the Boseglass-superfluid transition, but here occurs in a drivendissipative system. The existence of a state with longrange coherence, which is robust against disorder, may be important for applications of polariton condensation in areas such as analog simulation.
The synchronization transition we predict is absent for the Kuramoto model: it arises from the non-odd coupling between the phases, which gives a relevant nonlinear term in the continuum limit. That same form will appear for any system of coupled oscillators,θ k = ω k + K <l> f (θ l − θ k ), in which the coupling function, f (x), is neither purely even nor odd. That will generally be the case, so that our work implies that many other coupled oscillator systems can support synchronized states, notwithstanding disorder in their frequencies.
We acknowledge funding from the Irish Research Council under award GOIPG/2019/2824. Some calculations were performed on the Lonsdale cluster maintained by the Trinity Centre for High Performance Computing. This cluster was funded by grants from Science Foundation Ireland.