Dynamics of oscillator populations with disorder in the coupling phase shifts

We study populations of oscillators, all-to-all coupled by means of quenched disordered phase shifts. While there is no traditional synchronization transition with a nonvanishing Kuramoto order parameter, the system demonstrates a specific order as the coupling strength increases. This order is characterized by partial phase locking, which is put into evidence by the introduced novel correlation order parameter, which is shown to grow monotonically with the coupling strength, and via frequency entrainment by following concentration of the oscillators frequencies. Simulations with phase oscillators, Stuart–Landau oscillators, and chaotic Roessler oscillators demonstrate similar scaling of the correlation order parameter with the coupling and the system size and also similar behavior of the frequencies with maximal entrainment (at which the standard deviation of the frequencies is reduced by a factor close to four) at some finite coupling.


I. INTRODUCTION
Synchronization phenomena in populations of oscillators have been widely studied since the pioneering works of Winfree [1] and Kuramoto [2].The Kuramoto model and its modifications, based on the phase dynamics, serve as paradigmatic models of synchronization (with a possibility, in many situations, of a fully analytical treatment).The relevance of the findings within phase dynamics is confirmed by studies of realistic models of limit-cycle oscillators, where the amplitude dynamics is also included.If limit cycles are close to a Hopf bifurcation point, a universal normal form description can be adopted, which is expressed as an ensemble of Stuart-Landau oscillators, which also exhibits the synchronization transition [3,4].Moreover, Kuramoto-type synchronization can be also observed in populations of chaotic oscillators [5].
Effects of disorder on the synchronization properties of large populations of coupled oscillators are a subject of intensive current research.Randomness can be restricted to properties of individual units; e.g., already in the original formulation of the Kuramoto model [2], a distribution of natural frequencies was assumed.Also, the interactions between the oscillators can hold randomness, e.g., in the form of a random (weighted or unweighted) network [6][7][8].In many cases, in the thermodynamic limit of an infinite number of units, a theory of synchronization transition can be constructed where the randomness is smeared out, and the averaged coupling, if attractive, leads to a synchronization transition [9].
A special class of problems appears if the random coupling is maximally frustrating, i.e., if it disappears after a naive averaging.To the best of our knowledge, the first model of such type has been suggested by Daido [10], who considered a population of oscillators with mutual couplings sampled from a Gaussian distribution with zero mean.Although such a model does not exhibit synchronization as a transition to a state with a non-vanishing order parameter, it demonstrates some entrainment properties.By exploiting an analogy with spin glasses [11], he found regimes with a slow relaxation to equilibrium.Furthermore, Daido found a volcano transition, where a distribution of random coupling inputs on oscillators changes its form from one with a maximum at zero to one with a maximum at a finite distance from zero.The Daido model and its modifications have been since intensively studied [12][13][14][15][16][17][18].
This paper focuses on another type of randomness in a population of oscillators, where couplings are characterized by random phase shifts [19].If the distribution of the random shifts is nonuniform, one obtains certain renormalized averaged interaction in the thermodynamic limit, which begets a standard synchronization transition [9].However, no standard synchronization transition occurs for a maximally frustrating disorder with a uniform distribution of the phase shifts in [0, 2π).As we will show in this paper, a certain ordered state nevertheless appears, and the present work aims to characterize it.
The paper is organized as follows.In Section II, we will formulate the three basic models that will be explored numerically: phase oscillators, Stuart-Landau limit cycle oscillators, and Roessler chaotic oscillators.In Section III, we introduce the tools by which the regimes will be characterized.In particular, we introduce a novel correlation order parameter and discuss its relation to the standard Kuramoto order parameter.Detailed numerical studies of populations of phase oscillators are presented in Section IV.Similar simulations for Stuart-Landau and Roessler oscillations reported in Section V show, on the one hand, that the effects found for phase oscillators are reproduced for these more realistic models, while, on the other hand, novel features related to the amplitude dynamics appear.We conclude in Section VI with a comparison to Daido's disorder model.

II. BASIC MODELS
In this section, we take three popular models of global coupling (Kuramoto-Sakaguchi model [20]; system of Stuart-Landau oscillators [3]; system of chaotic Roessler oscillators [5]), widely used for exploration of the collective dynamics and of the synchronization transition, and extend them to the case of disorder in the phase shifts.

A. Phase oscillators
The standard Kuramoto-Sakaguchi model of N oscillators with all-to-all couplings reads Here φ k are the phases of the oscillating units, and ω k are their natural frequencies.Parameter f is the coupling strength (for future compatibility, we do not normalize it by N ).Here, all mutual coupling terms are assumed to include the same phase shift α.We introduce a disorder in these phase shifts by assuming that phase shifts α jk are independent, identically distributed random numbers (cf.Ref. [19]): Furthermore, in this paper we assume the maximally possible disorder, where α jk are uniformly distributed in [0, 2π) (for non-uniform distributions, see theory in Ref. [9]).Furthermore, we mostly concentrate on the asymmetric case α jk ̸ = α kj .
The main physical motivation behind the disordered phase shifts is that such shifts can be associated with time delays in coupling.In particular, if the frequencies of all oscillators are close to a central frequency ω 0 , then the shifts can be expressed via the coupling time delays τ jk as α jk = ω 0 τ jk , where ω 0 is the central frequency.The basic assumption is that the time delays are smaller than the characteristic time of slow variations of the phases; see Ref. [21] for a precise formulation.Uniformity (or almost uniformity) of the distribution of the phase shifts corresponds to a broad distribution of the time delays (a situation where delays in coupling obey a gamma distribution has been treated in Ref. [22]).We mention here that in addition to time delays in coupling, other factors may influence the phase shifts.For example, transmitting a signal from the driving to the driven unit may include filters with random characteristics.
We mention here that following the pioneering work of Daido [10], one often considers an implementation of disorder in the coupling, where the randomness is not in the phase shifts but in the coupling constants, see, e.g., Refs.[12][13][14][15][16][17][18].We discuss a possible relation to these models in Section VI; here we only mention that the phenomenologies in the cases of the coupling strength disorder and the phase shift disorders are different.
Below, we consider a Gaussian distribution of natural frequencies ω k .Since one can rescale time in Eq. (2), we will set the variance of this distribution to one.Furthermore, to avoid potential finite-size effects due to realizations of the frequencies, we assign ω k not via random sampling, but via deterministic sampling using quantiles of the Gaussian distribution (see discussion on the differences between random and deterministic samplings in Ref. [23]).Below, we will also consider a version of Eq. ( 2) where all the oscillators have the same natural frequency (which is set to zero, ω k = 0), but are subject to independent Gaussian white noise forces: Here again, the intensity of noise is fixed via the time rescaling.
In the case of deterministic identical oscillators, one can rescale the coupling parameter to one, and the problem is reduced to where the only parameter is the size of the ensemble N .

B. Stuart-Landau oscillators
A Stuart-Landau oscillator (SLO) is a prototypic model of a self-sustained oscillator, which includes both the phase and the amplitude dynamics.We consider in this paper a network of N SLOs with random phase shifts in coupling Here a k are complex amplitudes of the units, and parameter µ defines the relaxation rate of the amplitudes (µ −1 is the characteristic relaxation time of the amplitude dynamics).In the limiting case of strongly stable amplitude, µ → ∞, the dynamics of the phases φ k = arg(a k ) follows Eq. ( 2); otherwise, the dynamics of the amplitudes influences that of the phases.

C. Chaotic Roessler-type oscillators
A chaotic Roessler oscillator [24] (RO) possesses a welldefined phase variable; and an ensemble of such sys-tems with a Kuramoto-Sakaguchi-type coupling Eq. ( 1) demonstrates a phase synchronization transition similar to the Kuramoto transition [5].Below we explore the effect of randomness in the phase shifts in the coupling.For this purpose, we slightly modify the original Roessler system to obtain equations close to the dynamics of the SLO: Here, a is the complex "amplitude" of the oscillations (related to "usual" variables of the Roessler system as a = x + iy) , and z is a variable, the dynamics of which ensures chaotic "saturation" of the unstable oscillations of a.We introduce a parameter ω responsible for the natural frequency as a common factor to avoid difficulties related to possible bifurcations "order-chaos" in dependence on such a parameter (cf. a similar approach in Ref. [25]).The phase variable for the Roessler-type oscillator can be defined as φ k = arg(a k ).The effective dynamics of the phase is diffusive due to chaotic variations of the amplitudes.We now introduce a global coupling in a population with different natural frequencies ω k like in ( 5)

III. CHARACTERIZATIONS OF THE COLLECTIVE DYNAMICS
There are two ways to characterize synchronization in populations of oscillators: To follow the phase locking or the frequency entrainment.In many situations, these aspects have significant overlap, but there are also cases where phase locking is accompanied by frequency antientrainment [26,27].We will see that the two characterizations deliver complementary characterizations for the abovementioned ensembles.

A. Correlation order parameter for phase locking
Phase locking in a population of oscillators without disorder is typically characterized by the complex Kuramoto order parameter (complex KOP) At the synchronization transition, the absolute value of Z starts to grow, and large values of |Z| indicate that the phases are concentrated in a small range (compared to 2π).As we will see below, this KOP is nearly vanishing (up to finite-size fluctuations) in the presence of disorder.Therefore, we suggest to characterize relations between the phases with a complex phase correlation order parameter (complex COP) defined for the systems Eqs. ( 2),( 3),(4) as We notice that the imaginary part of Q k is, up to a factor, the force acting on the oscillator k.The real part of Q k will be more important.This quantity characterizes the closeness (on average) of the phase φ k to the value φ j + α jk , which would be the locked state if only one coupling j → k would be present (but this value is, of course, frustrated in the presence of all couplings).
Let us first demonstrate that complex COP and complex KOP are directly related in the case without disorder α jk = 0, i.e., for the standard Kuramoto setup: For large N , one thus has |Z| 2 ≈ Q (note that in the ordered case, Q is always real).Below, in the applications of the complex COP to the disordered populations, we will perform additional averaging of Q over time (for a particular sample of phase shifts), and over realizations of disorder, and separate real and imaginary parts of the resulting complex quantity as We will refer to C as the real correlation order parameter (real COP).Correspondingly, in a slight discrepancy with a usual definition of the (real) Kuarmoto order parameter, we will define it as |Z| 2 , averaged over time and over the realizations of disorder: The COP Q can be directly applied to populations of SLOs (5) and ROs (7), provided the phases exp[iφ k ] = a k /|a k | are used as the observables.One can also define the corresponding correlation order parameters that include full amplitudes where a j are complex amplitudes either of SLOs or of ROs.Similar to (10), Q (a) = C (a) + iD (a) .

B. Frequency entrainment
The definition of the average observed frequency for each oscillator in the population is straightforward (of course, the average includes only time average, not averaging over the realizations of disorder): This definition will be applied to all types of oscillators (phase oscillators, SLOs, and ROs).In the Kuramoto system without disorder, the spread of frequencies decreases beyond the synchronization transition.Moreover, in the deterministic case, a portion (growing with the coupling strength) of oscillators is perfectly entrained so that their frequencies coincide.
We notice here that the difference between the observed frequency Ω k (13) and the natural frequency ω k in (2) can be expressed as the average of the imaginary part of the local complex COP ( We also note that the entrainment of frequencies is a relevant characterization if the intrinsic frequencies ω k are different; for the noisy identical oscillators without the disorder, the frequencies are always equal.

IV. PHASE OSCILLATORS
In this section, we apply the introduced above characterizations of synchrony in a population of oscillators to the phase oscillators models Eqs. ( 2),(3),(4).

A. Phase oscillators with a Gaussian distribution of frequencies
In exploring model (2) for different system sizes N , we select the frequencies according to a Gaussian distribution with unit variance deterministically as corresponding quantiles.This is done to avoid additional diversity due to the random sampling of frequencies.The only parameters of the model are the coupling strength f and the system size N .The phase shifts in each run are chosen randomly according to a uniform distribution 0 ≤ α jk < 2π, an averaging is performed over time and over different samples of phase shifts.
Figure 1 presents the characterization in terms of the KOP and of the COP.The main finding is that while the value of the KOP R defined in Eq. ( 11) exhibits finite size fluctuations ∼ N −1 , which are practically independent on the coupling strength f (panel (c)), the value of the real COP C defined in Eq. ( 10) exhibits a monotonic growth with a saturation at large f (panel (a)).At the same time, the imaginary COP D fluctuates close to zero (panel (c)).
In Fig. 1(b) we also demonstrate, that the COP C fulfills a scaling relation Remarkably, the COP demonstrates a continuous transition without any particular threshold value of the coupling.
Next, in Fig. 2, we explore the properties of the frequency entrainment.We show in panel (a) the dependencies of frequencies Ω k according to Eq. ( 13) for one realization of disorder in dependence on the coupling strength f , and in panel (b) the standard deviations of the frequencies from the mean value, averaged for several realizations of disorder.One can see that at small coupling strengths f , the dispersion of frequencies decreases (from the value 1, defined according to the distribution of natural frequencies), but after reaching a minimum, it starts to increase again, nearly linear in f .We will discuss this feature and explain the value of the slope of the dashed line in more detail in Section IV C below.

B. Noisy phase oscillators
Here, we apply the same analysis as in Section IV A above to noisy oscillators described by Eqs.(3).First, we show in Fig. 3 the results for the order parameters characterizing phase coherence.One can see that there are no qualitative differences with the case of deterministic oscillators (Fig. 1).However, there are differences in the behavior of frequencies.Since the oscillators are identical, the observed frequencies Ω k in the absence of coupling coincide.However, with the onset of the coupling, such differences do appear, and, as Fig. 4 shows, the standard deviation grows linearly with f .This amplification corresponds to the increase of the dispersion of the frequencies for non-identical oscillators Eq. ( 2) at large coupling strengths f (Fig. 2(b)).

C. Identical phase oscillators
The results of Sections IV A, IV B above show that at large coupling strengths f , a universal regime appears, at which the real COP saturates at a value C sat ≈ 0.8N −1/2 , and the dispersion of the frequencies grows proportional to f .In this state, one can neglect differences between the oscillators (either due to their different natural frequencies or due to independent noise terms) and study a population of identical deterministic oscillators as defined in Eq. ( 4).Because now the time is rescaled according to the coupling strength, a constant dispersion of the observed frequencies will explain linear growth of this dispersion with f in setups where the coupling strength is not rescaled, i.e., in Figs.2(b),4.
The only parameter of Eqs. ( 4) is the population size N , what makes a statistical evaluation of different characteristics easier.First, we would like to characterize the dynamical regimes.Computation of the maximal Lyapunov exponent for the deterministic system (4) allows us to distinguish chaos and regular states.Because the system ( 4) is invariant with respect to the shift of all phases, one Lyapunov exponent is always zero.To distinguish whether a state with vanishing maximal Lyapunov exponent is a synchronous one where all the frequencies coincide or a quasiperiodic state where the frequencies are different, we calculated in each run the spread of the observed frequencies (13) as ∆Ω = Ω max − Ω min .Regimes where this spread nearly vanishes (we used a threshold 10 −4 as a practical criterion) were identified as steady states (in some reference frame rotating with a common frequency of the oscillators).Otherwise, the states with nearly vanishing maximal Lyapunov exponent and ∆Ω > 10 −4 were identified as quasiperiodic regimes.The probabilities of observing three possible regimes depend on the ensemble size as depicted in Fig. 5.For N = 2, for all phase shifts, the oscillators synchronize (as the theory of coupled oscillators predicts, because there is no frustration for two oscillators).At the same time, for larger N , quasiperiodic and chaotic regimes become possible.Finally, for large N ≳ 100, chaotic regimes dominate.
Next, we explored the statistics of the observed frequencies Eq. ( 13) at large N , and found that the distribution of these frequencies (when averaged over many realizations of disorder) is Gaussian with high accuracy.Furthermore, the variance is practically the same in a large range of N .This observation explains the linear behavior of the standard deviations of the frequencies in Figs.2(b) and 4, and also the fact that these standard deviations for different N nearly coincide at a value ≈ 0.274.This gives a prediction Std(Ω) ∼ 0.274f for large f in the presence of a distribution of natural frequencies or in the presence of noise; this prediction is shown with dashed lines in Figs.2(b) and 4.
Finally, we studied the statistics of the complex COP Q k defined in Eq. ( 8).We have found that the values Re(Q k ) and Im(Q k ) are practically uncorrelated.The average of Im(Q k ) is nearly zero, while the average of Re(Q k ) gives the averaged COP C, as discussed above.The variances of Re(Q k ) and Im(Q k ) are found to be ∼ N −1 .Remarkably, the distribution of Re(Q k ) is rather close to a Gaussian one, with small skewness values and excess kurtosis values.The distribution of Im(Q k ) is nearly symmetric with a small skewness, but the excess kurtosis values of ≈ 3 indicate deviations from a Gaussian distribution.

V. STUART-LANDAU AND ROESSLER-TYPE OSCILLATORS
This section aims to show that the essential effects described above for the phase oscillators are also observed for SLOs and ROs.We do not repeat all the numerical setups but focus on ensembles with a distribution of natural frequencies without noise.

A. Stuart-Landau oscillators
Stuart-Landau oscillators have a parameter µ that governs the stability of the limit cycle.For µ → ∞, a reduc-tion to the phase dynamics Eq. ( 2) becomes exact.Thus, we focus mainly on additional effects appearing at small µ.We start with presenting in Fig. 6 the evolution of the frequencies for an ensemble of SLOs (5), where natural frequencies are sampled according to a Gaussian distribution with variance 1.One can see that the evolution of standard deviation is very similar to that for phase oscillators (cf.Fig. 2): entrainment at small values of coupling and growth of the standard deviation for strong coupling.The dependence of these effects on the parameter µ is only quantitative: for small values of µ, the minimal standard deviation is achieved at smaller values of f (see in Fig. 6 the curves with an open square marker, corresponding to µ = 0.1).As discussed above, for the COP, we now have two variants: one is based on the transformation to the phases and thus is mostly close to the COP for the phase oscillators (Eq.( 8)); another one is based on the cross-correlations of the amplitudes (Eq.( 12)).Both these COPs are depicted in Fig. 7, and they demonstrate very different behaviors.While the phase-COP is practically µ-independent, the amplitude-COP grows strongly with the coupling if µ is small.We explain this by noticing that the linear coupling of the SLOs affects the stability of the equilibrium a k = 0 in a population of SLOs.Indeed, Eqs. (5) for small amplitudes a k can be written as a linear system and stability is affected by the random coupling matrix e −iα jk .According to the circular law [28], for large N , the eigenvalues of this matrix are uniformly distributed in a disk with radius √ N around the origin.Thus, the largest real part of the eigenvalues in Eq. ( 15) can be estimated as μ ≈ µ + f √ N .If we assume that this instability saturates at amplitude |a| 2 ≈ μ/µ, we conclude that the average amplitude of oscillations, for large couplings, grows ∼ f √ N /µ.This linear growth corresponds to what is observed in Fig. 7, and explains why the amplitude-COP is not a proper quantity to characterize phase coherence in the ensemble for small values of µ.

B. Roessler oscillators
Here, we report on the effects of coupling with disorder in the phase shifts on the Roessler system (7).The natural frequencies ω k have been sampled from Gaussian distribution with mean value 1 (because the RO is not rotationally invariant, one cannot set the mean frequency to zero by a transformation to a rotating reference frame, as it is possible for the phase oscillators and the SLOs) and standard deviation 0.1.Other parameters were set to d = 0.09, b = 0.4, c = 8.5.We also note that because the disordered coupling enhances instability of the oscillations, for large enough values of f , no finite regimes in system (7) have been observed because a trajectory escaped to infinity.Thus, we report on the regimes for relatively small values of coupling strength f only.
The evolution of the standard deviations of the observed frequencies with f is reported in Fig. 8.One can see that in the explored range of system sizes, this standard deviation is a universal function of rescaled coupling f √ N .Frequency entrainment is maximal at f √ N ≈ 0.14; however, the overall decrease of the standard deviation is less pronounced than for the phase oscillators and the SLOs (cf.Fig. 2(b) and Fig. 6).
Calculations of the COP, depicted in Fig. 9, revealed the same scaling in dependence on the f, N as in relation (14).The main difference between the corresponding results for the phase oscillators and the SLOs is that the imaginary part D also deviates from zero.We attribute this to the fact that the RO Eq. ( 7), due to the existence The difference in the amplitude of the values of Q and Q (a) has the same explanation as in the case of the SLOs above: disordered coupling contributes to the instability of oscillations of the variable a, thus increasing the amplitude of the chaotic attractor, which contributes to the growth of Q (a) . •

VI. CONCLUSION
We first summarize the results of our study.We have demonstrated that globally coupled oscillators with a maximally frustrating disorder in the coupling phase shifts demonstrate a transition to a partially ordered state.This state cannot be identified by the standard Kuramoto order parameter, which remains close to zero.We introduced a correlation order parameter, the calcula-tion of which requires knowledge of the phase shifts, and demonstrated that it continuously grows with coupling achieving a saturation level ∼ N −1/2 beyond a characteristic coupling strength ∼ N −1/2 (expression ( 14)).
The transition to order can also be characterized via frequency entrainment; here, one does not need to know the phase shifts, and this characterization can be easily implemented in experimental setups.Coupling leads to a concentration of frequencies, which can be straightforwardly quantified by a reduced standard deviation of their distribution.Remarkably, there is an optimal coupling strength at which the concentration is maximal; beyond this strength, the frequencies start to spread again.The disorder-induced spread of the frequencies in a population of identical oscillators explains the latter effect.
We have demonstrated that the ordering behavior described above holds for oscillators with a distribution of natural frequencies and noisy identical oscillators.Furthermore, we showed the effect not only for the phase oscillators but also for the Stuart-Landau oscillators and for the coupled chaotic Roessler systems.We focused on the case where the phase shifts are not symmetric (i.e., α jk ̸ = α kj ).Selected simulations with symmetric random phase shifts (i.e., α jk = α kj ) demonstrated only minimal quantitative difference to the presented features.
We have explored populations of oscillators with sizes up to several hundred; the scaling properties are wellestablished in this range.However, probably not all features can be downscaled to small ensembles.The reason is that for small populations, regular quasiperiodic (and periodic for identical deterministic units) become dominant, as the calculations of the largest Lyapunov exponent show.However, for N ≳ 100, chaotic regimes dominate.
Next, we compare our findings with other results on the disordered populations of coupled oscillators from the recent literature.We have found only paper [19] where a Kuramoto-type model with random phase shifts has been explored.Adopting the replica method for noisy phase oscillators with a Gaussian distribution of natural frequencies, the authors predicted for a maximally frustrating disordered case (phase shifts uniformly distributed in [0, 2π)) the appearance of a glass synchronization state if f N 1/2 exceeds a certain level, which nontrivially depends on the noise intensity.It is unclear if this "glass synchronization" state corresponds to a saturated ordered regime observed in our study at large couplings (no numerics is presented in Ref. [19]); such a comparison deserves further investigation.
In this respect, we also mention that we have not found any indication of the extremely slow dynamics, typically associated with a "glassy" state.Such a slowing has been observed in the Kuramoto model with another type of disorder, not in the phase shifts, but in the coupling constants (which are sampled from a Gaussian distribution with zero mean), see Refs.[10,[12][13][14][15][16][17][18].In the latter model also a sharp "volcano transition" has been observed; we have not found this transition in the random phase shifts model.This comparison shows that different types of disorder apparently lead to different types of the dynamics.Thus, it is appealing to look at combinational models, where both types of the disorder (in the phase shifts and the coupling constants) are present; this work is in progress.

FIG. 1 .
FIG. 1. Panel (a): real COP C (see Eqs. (10),(8)) for different system sizes in dependence on f .Panel (b): the same data as in panel (a), but in rescaled coordinates.Panel (c): imaginary COP D (see Eqs. (10),(8); scattered open markers in the bottom part of the panel) and R (see Eq. (11), filled markers at the top part of the panel; the values for different f nearly overlap) for different N .The dashed line shows the law ∼ N −1 for the values of R.

FIG. 2 .
FIG. 2. Panel (a): Evolution of the frequencies Ω k (see Eq. (13)) with coupling f for one realization of disorder (N = 100).Panel (b): standard deviation of the frequencies (averaged over many realizations) vs f for different N .The straight dashed line has a slope of 0.274 (see discussion in text).

FIG. 7 .
FIG. 7. Ensembles of SL oscillators.Panel (a): read COP C according to Eq. (8) vs. f , in rescaled coordinates.Panel (b): real COP C (a) according to Eq. (12) vs f , in rescaled coordinates.The values of N, µ are marked like in Fig.6.

FIG. 8 .
FIG. 8. Standard deviations frequencies Ω k vs f √ N for one realization of disorder and different N , for Roessler oscillators