Dynamics of oscillator populations with disorder in the coupling phase shifts

The standard Kuramoto-Sakaguchi model of N oscillators with all-to-all couplings reads

$$\begin{equation} \dot\varphi_k = \omega_k+f\sum_{j\neq k} \text{sin}\left(\varphi_j-\varphi_k-\alpha\right),\quad k = 1,\ldots,N\;. \end{equation} \tag{ 1 }$$

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 [21]):

$$\begin{equation} \dot\varphi_k = \omega_k+f\sum_{j\neq k} \text{sin}\left(\varphi_j-\varphi_k-\alpha_{jk}\right)\;. \end{equation} \tag{ 2 }$$

Furthermore, in this paper we assume the maximally possible disorder, where α_jk are uniformly distributed in $[0,2\pi)$ (for non-uniform distributions, see theory in [9]). Furthermore, we mostly concentrate on the asymmetric case $\alpha_{jk}\neq \alpha_{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 ω₀, then the shifts can be expressed via the coupling time delays τ_jk as $\alpha_{jk} = \omega_0\tau_{jk}$, where ω₀ 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 [23] 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 [24]). 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 [11], 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. [13–19]. We discuss a possible relation to these models in section 6; 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 equation (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 [25]). Below, we will also consider a version of equation (2) where all the oscillators have the same natural frequency (which is set to zero, $\omega_k = 0$), but are subject to independent Gaussian white noise forces:

$$\begin{equation} \begin{gathered} \dot\varphi_k = \eta_k\left(t\right)+f\sum_{j\neq k} \text{sin}\left(\varphi_j-\varphi_k-\alpha_{jk}\right)\;,\\ \left\langle \eta_k \right\rangle = 0,\;\left\langle \eta_k\left(t\right)\eta_l\left(t^{^{\prime}}\right) \right\rangle = 2\delta_{kl}\delta\left(t-t^{^{\prime}}\right)\;. \end{gathered} \end{equation} \tag{ 3 }$$

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

$$\begin{equation} \dot\varphi_k = \sum_{j\neq k} \text{sin}\left(\varphi_j-\varphi_k-\alpha_{jk}\right)\;, \end{equation} \tag{ 4 }$$

where the only parameter is the size of the ensemble N.

A 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

$$\begin{equation} \dot{a}_k = i\omega_k a_k +\mu a_k\left(1-|a_k|^2\right)+f\sum_{j\neq k} e^{-i\alpha_{jk}} a_j\;. \end{equation} \tag{ 5 }$$

Here a_k are complex amplitudes of the units, and parameter µ defines the relaxation rate of the amplitudes (µ⁻¹ is the characteristic relaxation time of the amplitude dynamics). In the limiting case of strongly stable amplitude, $\mu\to\infty$, the dynamics of the phases $\varphi_k = \text{arg}(a_k)$ follows equation (2); otherwise, the dynamics of the amplitudes influences that of the phases.

A chaotic Roessler oscillator [26] (RO) possesses a well-defined phase variable; and an ensemble of such systems with a Kuramoto–Sakaguchi-type coupling equation (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:

$$\begin{equation} \begin{aligned} \dot a& = \omega\left(i a+d a-z\right)\;, \\ \dot z& = \omega\left(b+z\left(\text{Re}\left(a\right)-c\right)\right)\;. \end{aligned} \end{equation} \tag{ 6 }$$

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 [27]). The phase variable for the Roessler-type oscillator can be defined as $\varphi_k = \text{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)

$$\begin{equation} \begin{aligned} \dot a_k& = \omega_k\left(i a_k+d a_k-z_k\right) +f \sum_{j\neq k} e^{-i\alpha_{jk}} a_j\;, \\ \dot z_k& = \omega_k\left(b+z_k\left(\text{Re}\left(a_k\right)-c\right)\right)\;. \end{aligned} \end{equation} \tag{ 7 }$$

Phase locking in a population of oscillators without disorder is typically characterized by the complex Kuramoto order parameter (complex KOP)

$$\begin{equation*} Z\left(t\right) = \frac{1}{N}\sum_{k = 1}^N e^{i\varphi_k\left(t\right)}\;. \end{equation*}$$

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 equations (2)–(4) as (cf [30, 31])

$$\begin{equation} \begin{gathered} Q\left(t\right) = \frac{1}{N}\sum_k Q_k\left(t\right)\;,\\ Q_k\left(t\right) = \frac{1}{N-1}\sum_{j\neq k} e^{i\left(\varphi_j-\varphi_k-\alpha_{jk}\right)}\;. \end{gathered} \end{equation} \tag{ 8 }$$

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 $\varphi_j+\alpha_{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 $\alpha_{jk} = 0$, i.e. for the standard Kuramoto setup:

$$\begin{equation} \begin{gathered} |Z|^2 = \frac{1}{N}\sum_{k = 1}^N e^{-i\varphi_k\left(t\right)}\frac{1}{N}\sum_{j = 1}^N e^{i\varphi_j\left(t\right)} = \frac{1}{N^2}\sum_{j,k} e^{i\left(\varphi_j-\varphi_k\right)} = \frac{\left(N-1\right)Q+1}{N}\;. \end{gathered} \end{equation} \tag{ 9 }$$

For large N, one thus has $|Z|^2\approx 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

$$\begin{equation} \left\langle Q \right\rangle = C+iD\;. \end{equation} \tag{ 10 }$$

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) Kuramoto order parameter, we will define it as $|Z|^2$, averaged over time and over the realizations of disorder:

$$\begin{equation} R = \left\langle |Z|^2 \right\rangle\;. \end{equation} \tag{ 11 }$$

The COP Q can be directly applied to populations of SLOs (5) and ROs (7), provided the phases $\text{exp}[i\varphi_k] = a_k/|a_k|$ are used as the observables. One can also define the corresponding correlation order parameters that include full amplitudes

$$\begin{equation} \begin{gathered} Q^{\left(a\right)}\left(t\right) = \frac{1}{N}\sum_k Q_k^{\left(a\right)}\left(t\right)\;,\\ Q_k^{\left(a\right)}\left(t\right) = \frac{1}{N-1}\sum_{j\neq k} a_j a_k^* e^{-i\alpha_{jk}}\;, \end{gathered} \end{equation} \tag{ 12 }$$

where a_j are complex amplitudes either of SLOs or of ROs. Similar to (10), $\left\langle Q^{(a)} \right\rangle = C^{(a)}+iD^{(a)}$.

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):

$$\begin{equation} \Omega_k = \left\langle \dot\varphi_k \right\rangle\;. \end{equation} \tag{ 13 }$$

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 $\Omega_k$ (13) and the natural frequency ω_k in (2) can be expressed as the average of the imaginary part of the local complex COP $(\Omega_k-\omega_k)\propto \left\langle \text{Im} (Q_k) \right\rangle$. 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.

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\unicode{x2A7D} \alpha_{jk}\lt2\pi$, 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 equation (11) exhibits finite size fluctuations ${\sim}N^{-1}$, which are practically independent on the coupling strength f (panel (c)), the value of the real COP C defined in equation (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)).

Zoom In Zoom Out Reset image size

In figure 1(b) we also demonstrate, that the COP C fulfills a scaling relation

$$\begin{equation} C\left(f,N\right) = N^{-1/2}c\left(f N^{1/2}\right)\;. \end{equation} \tag{ 14 }$$

Remarkably, the COP demonstrates a continuous transition without any particular threshold value of the coupling.

Next, in figure 2, we explore the properties of the frequency entrainment. We show in panel (a) the dependencies of frequencies $\Omega_k$ according to equation (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 4.3 below.

Zoom In Zoom Out Reset image size

Here, we apply the same analysis as in section 4.1 above to noisy oscillators described by equation (3). First, we show in figure 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 (figure 1). However, there are differences in the behavior of the frequencies. Since the oscillators are identical, the observed frequencies $\Omega_k$ in the absence of coupling coincide. However, with the onset of the coupling, such differences do appear, and, as figure 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 equation (2) at large coupling strengths f (figure 2(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results of sections 4.1 and 4.2 above show that at large coupling strengths f, a universal regime appears, at which the real COP saturates at a value $C_{sat}\approx 0.8 N^{-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 equation (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 figures 2(b) and 4.

The only parameter of equation (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 $\Delta\Omega = \Omega_\textrm{max}-\Omega_\textrm{min}$. Regimes where this spread nearly vanishes (we used a threshold 10⁻⁴ 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 $\Delta\Omega\gt10^{-4}$ were identified as quasiperiodic regimes. The probabilities of observing three possible regimes depend on the ensemble size as depicted in figure 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\gtrsim 100$, chaotic regimes dominate.

Zoom In Zoom Out Reset image size

Next, we explored the statistics of the observed frequencies equation (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 figures 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 $\text{Std}(\Omega)\sim 0.274 f$ 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 figures 2(b) and 4.

Finally, we studied the statistics of the complex COP Q_k defined in equation (8). We have found that the values $\text{Re}(Q_k)$ and $\text{Im}(Q_k)$ are practically uncorrelated. The average of $\text{Im}(Q_k)$ is nearly zero, while the average of $\text{Re}(Q_k)$ gives the averaged COP C, as discussed above. The variances of $\text{Re}(Q_k)$ and $\text{Im}(Q_k)$ are found to be ${\sim}N^{-1}$. Remarkably, the distribution of $\text{Re}(Q_k)$ is rather close to a Gaussian one, with small skewness values and excess kurtosis values. The distribution of $\text{Im}(Q_k)$ is nearly symmetric with a small skewness, but the excess kurtosis values of ${\approx}3$ indicate deviations from a Gaussian distribution.

SLOs have a parameter µ that governs the stability of the limit cycle. For $\mu\to\infty$, a reduction to the phase dynamics equation (2) becomes exact. Thus, we focus mainly on additional effects appearing at small µ.

We start with presenting in figure 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 figure 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 figure 6 the curves with an open square marker, corresponding to µ = 0.1).

Zoom In Zoom Out Reset image size

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 (equation (8)); another one is based on the cross-correlations of the amplitudes (equation (12)). Both these COPs are depicted in figure 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, equation (5) for small amplitudes a_k can be written as a linear system

$$\begin{equation} \dot{a}_k = i\omega_k a_k +\mu a_k+f\sum_{j\neq k} e^{-i\alpha_{jk}} a_j\;, \end{equation} \tag{ 15 }$$

and stability is affected by the random coupling matrix $e^{-i\alpha_{jk}}$. According to the circular law [32], for large N, the eigenvalues of this matrix are uniformly distributed in a disk with radius $\sqrt{N}$ around the origin. Thus, the largest real part of the eigenvalues in equation (15) can be estimated as $\tilde\mu\approx \mu+f\sqrt{N}$. If we assume that this instability saturates at amplitude $|a|^2\approx \tilde\mu/\mu$, we conclude that the average amplitude of oscillations, for large couplings, grows ${\sim}f\sqrt{N}/\mu$. This linear growth corresponds to what is observed in figure 7, and explains why the amplitude-COP is not a proper quantity to characterize phase coherence in the ensemble for small values of µ.

Zoom In Zoom Out Reset image size

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 figure 8. One can see that in the explored range of system sizes, this standard deviation is a universal function of rescaled coupling $f\sqrt{N}$. Frequency entrainment is maximal at $f\sqrt{N}\approx 0.14$; however, the overall decrease of the standard deviation is less pronounced than for the phase oscillators and the SLOs (cf figures 2(b) and 6).

Zoom In Zoom Out Reset image size

Calculations of the COP, depicted in figure 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 equation (7), due to the existence of the variable z, is non-invariant toward rotations of the phase $a\to a e^{i\theta}$.

Zoom In Zoom Out Reset image size

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)}$.