Perfect synchronization in networks of phase-frustrated oscillators

Synchronization captures the emergence of collective behavior in complex systems [1–3], ranging from neuronal dynamics [4] to animal behavior [5] and technological networks [6]. In its classic formulation synchronization is driven by the coupling between the oscillators, which drives them towards collective oscillations, overcoming the diversity in the intrinsic frequencies of each individual oscillator [7–11]. Hence synchronization is enhanced either by increasing the coupling strength between the oscillating units, or by a homogeneous frequency distribution among all oscillators. These strategies towards synchronization, however, fail in case the coupling between the oscillators induces phase lags [12–14], a common characteristic featured by many real systems, where the components take time to respond to their neighboring oscillators. Indeed, under phase-frustration, the system persistently avoids synchronization, even when the frequencies are homogeneous or under relatively strong coupling.

To overcome this challenge, we derive here the link between the network characteristics, the phase lags and the optimal frequency set, that allows the phase-frustrated system to reach perfect synchronization. This allows us a two-way prediction: for a given network and phase lags, we predict the optimal selection of natural frequencies that will ensure synchronization. Alternatively, given a set of natural frequencies, we show how to design the network that will lead the oscillators towards perfect synchrony. We find, numerically, that our predicted synchronization is quite robust, exhibiting a range of phase-locked solutions even under deviations from our predicted frequencies/networks, thus being insensitive to moderate levels of noise/perturbation. Counterintuitively, we find that synchronization is not necessarily enhanced by strengthening the coupling or by selecting homogeneous frequencies, rather it emerges from the complex interplay between the selected frequencies, the distribution of phase lags and the structure of the weighted underlying network.

Consider a system of N coupled oscillators, whose phases $\theta_i(t)$, $(i = 1,\ldots,N)$, are driven by the dynamic equations [7,12]

$$\begin{equation} \frac{\text{d}\theta_i}{\text{d}t} = \omega_i + \sum_{j=1}^{N} A_{ij}{F}(\theta_j - \theta_i-\alpha_{ij}), \end{equation} \tag{ 1 }$$

where ω_i represents node i's natural frequency and A_ij is a weighted adjacency matrix with arbitrary degree and weight distributions. The functional form of the coupling is captured by $F(\theta_j - \theta_i - \alpha_{ij})$, a $2 \pi$ periodic function, with distributed phase lags $\alpha_{ij}$, which capture the response time of oscillator i to changes in its neighbor's phase $\theta_j$. Phase-frustrated models of the form (1) are frequently used to describe coupled systems, from Josephson junctions [15] to power supply networks [16] and mechanical rotors [17]. Choosing $F(\theta) = \sin(\theta - \alpha)$ (with $\alpha_{ij}=\alpha$ independent of i and j), eq. (1) converges to the Sakaguchi-Kuramoto model of phase frustrated oscillators [12–14,18]; setting $\alpha \rightarrow 0$ we arrive at the classic Kuramoto dynamics [7].

To quantify the level of synchronization in the system we use the Kuramoto order parameter

$$\begin{equation} r = \left|\frac{1}{N}\sum_{j=1}^{N} e^{i\theta_j}\right|, \end{equation} \tag{ 2 }$$

which approaches $r = 0$ in the disordered regime and r = 1 in the limit where all oscillators are in perfect synchrony, namely $\theta_1 = \theta_2 = \dots = \theta_N$. In the classic Kuramoto framework the level of synchronization is determined by the trade-off between the heterogeneity of the natural frequencies ω_i and the strength of the coupling, as determined by A_ij. Hence to achieve synchronous behavior one draws ω_i from a narrowly bounded distribution (e.g., normal distribution) or increases the weights of A_ij until reaching $r \rightarrow 1$. Such perfect synchronization, however, is unattainable in the presence of phase-lags even for extreme coupling strengths [19,20]. Hence we seek the optimal frequency sequence ${\boldsymbol \omega} = (\omega_1, \omega_2,\ldots.\,,\omega_N)^\top$ that will enable perfect synchronization r = 1 for phase-frustrated oscillators of the form (1).

To obtain ω, we analyze eq. (1) as it approaches synchronization [11,19–21], namely in the limit where ${|\theta_j - \theta_i| \rightarrow 0}$. In this limit, the coupling function in (1) can be approximated by $F(\theta_j - \theta_i - \alpha_{ij}) \approx F(-\alpha_{ij}) + F^{\prime}(-\alpha_{ij})(\theta_j - \theta_i)$, where $F^{\prime}(\alpha) = \left. \frac{\partial F}{\partial\theta} \right|_{\theta = 0}$. This allows us to write eq. (1) as

$$\begin{equation} \frac{\text{d}\theta_i}{\text{d}t}= \omega_i + d_{i}- \sum_{j=1}^{N}L_{ij}\theta_j, \end{equation} \tag{ 3 }$$

where

$$\begin{equation} {{d}_i} = \sum_{j=1}^{N}A_{ij} F(-\alpha_{ij}) \end{equation} \tag{ 4 }$$

and

$$\begin{equation} {L}_{ij}=\delta_{ij}\left(\sum_{q=1}^{N} A_{iq}F^{\prime}(-\alpha_{iq})\right)- A_{ij} F^{\prime}(-\alpha_{ij}), \end{equation} \tag{ 5 }$$

in which $\delta_{ij}$ is the Kronecker δ-function. The system will reach a synchronized state if, for some choice of the natural frequencies ω, eq. (3) reaches a stable solution in which all phases, $\theta_i(t)$, evolve according to some common frequency Ω, namely $\theta_i(t) = \phi_{i}+\Omega t$, where $\phi_{i}=\theta_i(0)$ is the phase of oscillator i. Transforming to an Ω-rotating frame, we have $\frac{\text{d}\theta_i}{\text{d}t}=0$, which in eq. (3) leads to

$$\begin{equation} {\boldsymbol {\phi}} = { L^\dagger}{({\bm{\omega}} + {\bm{d}})}, \end{equation} \tag{ 6 }$$

where ${\boldsymbol {\phi}}=(\phi_{1},\phi_{2},\ldots,\phi_{N})^\top$ represents the vector of all phases, ${\boldsymbol {d}}=(d_{1},d_{2},\ldots,d_{N})^\top$ and ${{ L}^\dagger}$ is the pseudo-inverse [22] of L in (5). The condition (6) captures frequency synchronization, a state in which all units oscillate at a common frequency, but with different phases $\phi_i$. Complete synchronization, however, requires also that all phases condense around a common value ϕ, which, by additional rotation of the system can be set to ${\boldsymbol {\phi}}=0$. With this gauging with arrive at

$$\begin{equation} {\bm{\omega}} = -{\bm{d}}+ C \mathbf{1}, \end{equation} \tag{ 7 }$$

where $\mathbf{1}=(1,\ldots,1)^\top$ and C is an arbitrary constant, a degree of freedom enabled due to the fact that ${L^{\dagger} \mathbf{1}=0}$ [22]. We use this degree of freedom to select $C= \langle{d}\rangle$, allowing us to write, explicitly, the designated frequency of the i-th oscillator as

$$\begin{equation} \omega_{i} = -\sum_{j=1}^{N} A_{ij}F(-\alpha_{ij})+\frac{1}{N} \sum_{i,j=1}^{N} A_{ij}F(-\alpha_{ij}), \end{equation} \tag{ 8 }$$

for which $\langle {\bm{\omega}}\rangle=\frac{1}{N} \sum_{i=1}^{N} \omega_{i}=0$, namely we gauge the mean frequency to be zero.

Equation (8) represents our key prediction, providing the optimal frequency set ω in a weighted network of heterogeneous phase-frustrated oscillators. It indicates that the optimal frequency set ω is determined by the interplay between the system's topology (A_ij), its dynamics (F) and the specific form of the distributed phase lags $(\alpha_{ij})$. For the unfrustrated Kuramoto model $(F(\theta)=\sin(\theta),\alpha=0)$ it predicts that the optimal frequency set is uniform, ${w_i = 0}$ for all i, reaffirming Kuramoto's classic prediction [7]. In the Sakaguchi-Kuramoto model $(\alpha_{ij}=\alpha)$ eq. (8) predicts that ω_i scales with node i's weighted degree $\omega_i \sim \sum_{j=1}^N A_{ij}$ up to an additive constant $C=\frac{F(-\alpha)}{N}\sum_{i,j=1}^N A_{ij}$. This implies that contrary to the Kuramoto model, where synchrony is a consequence of ω 's homogeneity, in the phase-frustrated case ω depends on A_ij's degree sequence, therefore it must follow A_ij's degree heterogeneity. For instance, if A_ij is scale-free, as often encountered in real networks [23,24], ω_i must also be drawn from a scale-free distribution. Hence, counterintuitively, (8) shows that perfect synchrony may arise from oscillator heterogeneity, namely from a scale-free sequence ω.

To test our prediction we constructed eq. (1) using a weighted scale-free network A_ij ($P(k)\sim k^{-\gamma}, \gamma=3$) of $N = 1,000$ interacting oscillators, whose phase lags were extracted from a uniform distribution $0.1 \le \alpha_{ij} \le \pi/2$, i.e. $\alpha_{ij} \sim \mathcal{U}(0.1,1.57)$. The weights of all existing links were also extracted from a uniform distribution ${A_{ij}\sim \mathcal {U}(0.1,1.5)}$. We then numerically solved eq. (1) and tested the level of synchronization r (2), for several choices of ω: homogeneous, in which all ω_i are identical; normal, in which ω_i are extracted from a normal distribution with mean 0 and variance 1, i.e. $\omega_{i}\sim \mathcal{N}(0,1)$; and uniform, where $\omega_{i}\sim \mathcal{U}(-2,2)$. For uniform ω we find that the system cannot synchronize with r being significantly smaller than unity (fig. 1(a), red). Synchronization becomes even lower for normal (blue), and slightly improved for homogeneous (green). Hence, as opposed to the unfrustrated Kuramoto dynamics, here a bounded frequency distribution cannot lead to synchronization. Our theory, however, predicts that synchronization can be obtained if we construct ω using the optimal frequency set (8). Indeed, we find that selecting our predicted ω, the system successfully reaches perfect synchronization, featuring r = 1, as predicted (black).

Zoom In Zoom Out Reset image size

As explained above, synchronization is often a consequence of two competing effects: the strength of the coupling, which forces the system into collective oscillations vs. the heterogeneity in ω, which drives the system away from synchronization. We now test these two effects systematically, by first, rescaling, and hence weakening/strengthening, the coupling between all oscillators, and then adding increasing levels of noise to the optimal frequency set predicted in (8).

Coupling strength: Often one wishes to force a system towards global synchrony by strengthening the level of coupling between the interacting oscillators, for instance, multiplying all A_ij terms (weights) by a factor of K > 1. For phase-frustrated systems of the form (1), however, such approach will not lead to global synchrony. Indeed, as fig. 1(b) indicates, for ω uniform (red), normal (blue) and homogeneous (green), the system consistently avoids global synchronization, despite increasing K. For the optimal frequency set (black), we obtain perfect synchronization for $K = 1$ (magenta dot), as predicted, yet increasing or decreasing K harms the level of synchronization since any change to A_ij, even increasing the strength of the coupling (K > 1), leads to consequent changes in the optimal frequency set (8). The important point is that for a rather broad range of K values, the system sustains relatively high levels of synchronization, allowing for a phase-locked solution, even if the selected frequency set is not precisely the optimal one predicted by (8). This represents the robustness of our solution, opening a wide window of phase-locked solutions in the vicinity of the optimal selection (8).

Frequency deviations: To test the sensitivity of the synchronization to deviations from the optimal frequency set (8) we add Gaussian noise to ω, setting ${\omega_i} \rightarrow {\omega_i} + {\delta{\omega_i}}$, where ${\delta{\omega_i}} \sim \mathcal {N}(0, \sigma{\omega_i})$, a normally distributed random variable with mean zero and variance ${\sigma^2{\omega_i}^2}$, representing multiplicative noise that is proportional to ω_i. Such deviation will reduce the level of synchronization to r < 1, resulting in synchronization loss, which can be quantified by $\rho = 1 - r\ (0 \le \rho \le 1)$. For small σ the deviation from synchronization is small, allowing us to approximate (2) up to second order as $r \approx 1 - \|{\bm{\phi}}\|^2/{2N}$ [11,19], therefore,

$$\begin{equation} \rho = 1 - r \sim \frac{1}{2} \rm {Var} (\phi), \end{equation} \tag{ 9 }$$

where Var(X) represents the variance of the random variable X. Using (6), we write

$$\begin{equation} {\phi_i} = { \sum_j L_{ij}^\dagger} {({\omega_j} + {d_j})} = \sum_j L_{ij}^\dagger\delta{\omega_j}, \end{equation} \tag{ 10 }$$

and hence, with $L_{ij}^\dagger$ being approximately independent of $\delta{\omega_j}$, we have

$$\begin{equation} {\rm{Var}} (\phi_{i})= \sum_j (L_{ij}^\dagger)^{2} {\rm{Var}}(\delta{\omega_{j}}), \end{equation} \tag{ 11 }$$

where we also used the fact that $\delta \omega_j$ are independent of each other. As a result we find that

$$\begin{equation} {\rm{Var}}(\phi_{i}) \approx C_{i} \sigma^2, \end{equation} \tag{ 12 }$$

where $C_{i}=\sum_{j}(L_{ij}^\dagger)^2 \omega_{j}^2$, and hence the overall variance of all $\phi_i$ satisfies ${\rm{Var}}(\phi) \propto \sigma^2$, where the proportion constant is a function of all pre-factors C_i. Omitting such factors that do not depend on the noise level σ, we arrive at the scaling relationship ${\rm Var(\phi)}\sim \sigma^2$, which in (9) provides

$$\begin{equation} \rho\sim\sigma^2, \end{equation} \tag{ 13 }$$

showing that synchronization loss is quadratically dependent on the noise level in the oscillator frequencies. This allows us to evaluate the decay in synchronization as we deviate from the optimal frequency set (8). For small σ we predict the scaling (13) and as σ increases ω continues to deviate from (8) until eventually synchronization is completely lost and $\rho\rightarrow 1$. This behavior is clearly observed in fig. 2, where we introduce increasing levels of noise to the optimal frequency set. We find that for small noise levels $\rho\sim\sigma^2$ (solid line) as predicted and for $\sigma \rightarrow 1$, a limit where ω is completely overridden by noise, synchronization is lost with $r \rightarrow 0$ and, consequently, ρ approaching unity.

Zoom In Zoom Out Reset image size

Second-order dynamics: Our formalism is also applicable beyond the limits of eq. (1). To show this we focus on second-order phase-frustrated Kuramoto dynamics, captured by

$$\begin{equation} \frac{\text{d}^2\theta_i}{\text{d}t^2} = P_i -\beta\frac{\text{d}\theta_i}{\text{d}t} + \sum_{j=1}^{\mathrm{N}}A_{ij}F(\theta_j - \theta_i-\alpha_{ij}), \end{equation} \tag{ 14 }$$

as frequently used to describe phase synchronization in power supply networks [6,25–28]. In (14) P_i represents the generated (P_i > 0) or consumed (P_i < 0) power, and β is the damping coefficient of the system components. To examine our formalism in an empirical setting we collected data from the northern European power grid [29], comprising $N = 236$ nodes and $E = 320$ links. We extracted the frustration terms from a uniform distribution, $\alpha_{ij} \sim \mathcal{U}(0,0.5)$. As before, we constructed the optimal frequency set ω and tested the level of synchronization r against varying levels of coupling K and noise σ, setting $\beta = 10$ and $\beta = 0.1$, to represent the limits of strong and weak damping. We find that also in the case of second-order dynamics, the empirical power grid network reaches global synchronization $(r = 1)$ for $K = 1$, as predicted. It exhibits a phase-locked solution, $r \lesssim 1$, at a range of K values around $K = 1$ (fig. 3(a), (c)). The loss of synchronization scales as $\rho \sim \sigma^2$, for small deviations from the optimal ω, with complete synchronization loss $(\rho \rightarrow 1)$ at $\sigma \approx 1$, the point where noise levels become comparable to the frequencies themselves (fig. 3(b), (d)).

Zoom In Zoom Out Reset image size

Designing networks for synchronization: Our theory, up to this point, focused on the selection of ω that will enable (1) to reach synchronization, namely, we begin with a given weighted network A_ij and phases $\alpha_{ij}$, for which we seek the optimal set of natural frequencies ω_i. Often, however, we are confronted with the opposite challenge: given a set of oscillators with natural frequencies ω, can we design a weighted network A_ij with lags $\alpha_{ij}$ that will drive the oscillators toward synchronization? As indicated by eq. (8) this reverse challenge is not as well-defined, allowing a broad degree of freedom to select the network, its link weights and the matching phase lags, hence for a given set ω, one can construct many synchronizable networks. To examine this systematically we consider the case where the phase lags were all set to $\alpha_{ij} = \alpha = 0.1$, leaving us with the degree of freedom to construct A_ij. Extracting the natural frequencies from a uniform distribution $\omega_i \sim \mathcal{U}(2,30)$, we constructed a weighted network that satisfies (8), by setting its weighted degrees d_i (4) to conform with the condition (7). As predicted, we find that the designed network leads to perfect synchronization r = 1 (fig. 4(a)); perturbing this network results in gradual synchronization loss (fig. 4(b)).

Zoom In Zoom Out Reset image size

Dimension reduction analysis: To assess the stability of our observed synchronization we use a collective coordinate approach [21,30,31], to analyze the behavior of the synchronized and phase-locked solutions in the rotating frame. In this approach, the instantaneous phase of each phase-locked oscillator is approximated by

$$\begin{equation} \theta_i(t)=\chi\omega_i, \end{equation} \tag{ 15 }$$

where $\chi = \chi(t)$ is a time-dependent function. A stable phase-locked solution emerges in the network when the error

$$\begin{equation} \epsilon_{i}(\chi) = \dot{\chi}\omega_i-\omega_i-\sum_{j=1}^{N} A_{ij} F(\chi(w_j-w_i)-\alpha_{ij}), \end{equation} \tag{ 16 }$$

is minimized. Such minimization is enabled if $\epsilon_{i}({\chi})$ is orthogonal to the tangent subspace of the solution space of eq. (15), spanned by ${\frac{\partial \theta_i}{\partial \chi}}=w_i$ [30,31]. Projecting $\epsilon_{i}(\chi)$ onto this restricted subspace and using the orthogonality condition we obtain a reduced one-dimensional differential equation for $\chi(t)$

$$\begin{equation} \frac{\text{d}\mathbf{\chi}}{\text{d}t}= g(\chi), \end{equation} \tag{ 17 }$$

where

$$\begin{equation} g(\chi)=1+\frac{1}{\sigma^2}\sum_{i=1}^{N}\omega_i\sum_{j=1}^{N} A_{ij}F(\chi(w_j-w_i)-\alpha_{ij}) \end{equation} \tag{ 18 }$$

and $\sigma^2=\sum_{j=1}^{N}{\omega_j}^2$. Equation (17) reaches a stable fixed point if for some choice of ω, eq. (18) satisfies $g(\chi)=0$ and $g^\prime(\chi) = \frac{\text{d} {g}}{\text{d}{\chi}}< 0$. Under these conditions we have $\chi(t) = \chi$ independent of time, which in (15) predicts a time-independent $\theta_i$, representing a phase-locked system, where all nodes oscillate at a common frequency, with their relative phases constant in time. If such fixed-point occurs for $\chi \rightarrow 0$ the phases (15) approach a single value $\theta_i \rightarrow 0$, representing a convergence to prefect synchronization $r \rightarrow 1$. In fig. 5(a) we show $g(\chi)$ vs. χ taking our predicted optimal ω from (8). We find that a stable fixed point ($g(\chi) = 0$, $g^{\prime}(\chi) < 0$) occurs at $\chi = 0$, precisely the predicted global synchronization obtained under our optimal ω (8). For the other three selections of ω (uniform, normal, homogeneous), we observed that $g(\chi)$ never crosses zero, indicating that stable synchronization is indeed unattainable (fig. 5(a), red, blue, green, comparing to fig. 1(a)).

Zoom In Zoom Out Reset image size

To test the impact of changes to the coupling K, we show, in fig. 5(b), the behavior of $g(\chi)$ under the optimal ω, for a selection of K values. As explained above, perfect synchrony occurs when $g(\chi) = 0$ and $g^\prime(\chi) < 0$ at precisely $\chi = 0$, conditions that are only satisfied for $K = 1$ (black), the value for which our optimal ω was constructed. Especially important is the critical $K = 0.78$ (yellow), representing the coupling where (17) assumes a stable fixed-point for the first time. This point represents the onset of a phase-locked solution, where r begins to rise above zero. Beyond this point r consistently increases as the crossing point ($g(\chi) = 0$) approaches closer to $\chi = 0$, eventually reaching perfect synchronization at $K = 1$, for which $g(0) = 0$. Indeed, in fig. 1(b) we observe the onset of synchronization at precisely $K = 0.78$, a reassuring consistency with fig. 5(b). For larger values of $K\ (K = 1.2, 1.5)$ the system continues to exhibit a stable phase-locked solution, indicated by the slow decline in fig. 1(b) for K > 1 and by the consistent crossing of $g(\chi) = 0$ in fig. 5(b) (green, cyan). Note, that for $K = 1.5$ the fixed-point condition is satisfied for two values of χ, indicated by the two crossing points $g(\chi) = 0$, both featuring a negative slope. In our analysis, however, we only regard the first crossing point, as phase-locking is only captured by (15) in the limit of small χ.

The onset of synchronization: The optimal set ω is designed for perfect synchronization r = 1 at a given weight $K = K_{\rm Opt}$ (set to unity in our analyses up to this point). However, the onset of synchronization occurs at a critical $K_c < K_{\rm Opt}$, a point where r begins to rapidly ascend from the chaotic regime $r = 0$ ($K = 0.78$ in our previous example, figs. 1(b) and 5(b)). In the Supplementary Material Supplementarymaterial.pdf we use mean-field analysis [32–34] to analytically predict K_c under homogeneous phase lags $\alpha_{ij} = \alpha$ as

$$\begin{equation} K_c(\alpha)=\frac{2 K_{\rm Opt}^3 {\sin^{3} \alpha} {\langle k\rangle} \cos \alpha}{\pi (\Omega_c+b)^2 P\left(\frac{\Omega_c+b}{K_{\rm Opt}\sin(\alpha)}\right)}, \end{equation} \tag{ 19 }$$

where the group angular velocity $\Omega_c$ satisfies

$$\begin{eqnarray} &&\pi\left(\frac{\Omega_c +b}{a}\right)^2 P\left(\frac{\Omega_c +b}{a}\right)\tan \alpha=\int_{q_{min}}^{\infty} \frac{q^2 P(q)}{aq-b-\Omega_c} \text{d}q, \end{eqnarray} \tag{ 20 }$$

in which $a = K_{\rm Opt} \sin \alpha$, $b = K_{\rm Opt} \sin \alpha \frac{\sum_{i=1}^{N}q_i}{N}$, q_min is the minimum degree over all nodes and P(q) is the probability density function of node frequencies. Equations (19) and (20), our final prediction, allow us, for a given A_ij, α and natural frequencies ω, to express K_c, the critical point of transition, in which synchronization begins to emerge. Together with $K_{\rm Opt}$ of eq. (8) these two points fully characterize the states of the system: chaotic ($r = 0$ for K < K_c, phase-locked (r > 0) at $K \ge K_c$ and optimal $(r = 1)$ at $K = K_{\rm Opt}$.

To observe this we constructed a scale-free A_ij (${N = 5000}$, $\langle k \rangle = 30$, $P(k) \sim k^{-\gamma}$, $\gamma = 3$) with homogeneous phase lags α. We matched this network with the appropriate optimal frequency sets ω, such that perfect synchronization occurs at $K_{\rm Opt} = 1$. In fig. 6(a) we show r vs. K, for $\alpha = 0.1$ (blue) and $\alpha = 0.5$ (red) finding that indeed, in both cases, r = 1 at the optimal $K = K_{\rm Opt} = 1$. The crucial point is that synchronization begins to appear at significantly lower values of K, at the critical K_c, where r begins to sharply incline. We next used (19) to calculate K_c in both cases, confirming our analytical predictions, which perfectly match the observed criticality (vertical dashed lines). Repeating the same experiment, this time setting $K_{\rm Opt} = 0.5$ we observe further confirmation of our prediction (fig. 6(b)).

Zoom In Zoom Out Reset image size

Understanding the phenomena of synchronization in networks has crucial applications in fields ranging from neuronal networks to power supply. These systems are often described by highly heterogeneous weighted networks and exhibit distributed lag times, a combination that rarely succumbs to analytical treatment [35–38]. Our analysis here has shown how to analytically construct the appropriate frequency sequence to ensure perfect synchronization, relating the optimal frequency set to the weighted network structure and to the distributed lags $\alpha_{ij}$. It shows that the systems' weighted degree distribution plays an important role in determining the desired frequencies, where degree heterogeneity dictates a similar heterogeneity in the frequency set. This shows that synchronization can occur by introducing diversity in ω_i, rather than by increasing its homogeneity. Hence cooperative phenomena may emerge even in the presence of microscopic diversity, a consequence of the phase lags, that is absent in the classical Kuramoto framework.

The authors would like to thank Syamal Dana for interesting comments and suggestions. PK acknowledges support from DST, India under the DST-INSPIRE scheme (Code: IF140880). CH is supported by the CHE/PBC, Israel. This work was supported by NSF-CRISP Award No. 1735505.