Frustrated Synchronization of the Kuramoto Model on Complex Networks

1. Introduction

Synchronization phenomena are abundant in nature ranging from brain science [1] to power grids [2]. Describing them via toy models has a long history, for recent reviews see [3,4]. One of the simplest and well known is introduced by Kuramoto [5]. In its original form it is a fully coupled system defined on a full graph. The synchronization transitions of the locally coupled versions were studied on finite dimensional lattices [6]. In homogeneous systems the Kuramoto model exhibits $d_c=d_l=4$, that means real phase transition can happen above $d\ge d_l=4$ and it is mean-field type. Entrainment transition of frequency variables can be non-mean-field like for $2<d<4$.

On heterogeneous, random graphs the phase transition remains mean-field like [7] according to the annealed heterogeneous mean-field approximation, however care must be taken even in dense-network systems, particularly in the disordered phase [8]. By studying the dynamical scaling of Kuramoto on Erdős Rényi graphs non-monotonic corrections to scaling hindered to see clearly mean-field criticality even on very large systems [9]. Note the disorder here is different from other ‘disordered’ Kuramoto models studied in the literature, where the disorder arises from independent random positive and negative couplings, which add frustration to the system and can lead to glassy dynamics [10].

In lower dimensional, simplicial complex model of manifolds complexes [11] and on hierarchical modular networks [12,13] a so-called frustrated synchronization transition was reported for spectral dimension $d_s<4$. Ref [14] showed that synchronized phase can only be thermodynamically stable for spectral dimensions above four and that phase entrainment of the oscillators can only be found for spectral dimensions greater than two. Very recently Ref. [15] studied one-dimensional long-range random ring networks, where any two nodes on the network are connected with a probability proportional to a power law of the distance between the nodes and confirmed the results of [14] for the Kuramoto model.

On a large, weighted human connectome network, containing 804 092 nodes, the topological dimension is $d<4$, a real synchronization phase transition is not possible in the thermodynamic limit, still a transition between partially synchronized and desynchronized states could be found, with non-universal coupling dependent dynamical scaling [16]. On the other hand on the infinite dimensional “2dll” random long link model of lateral size $L=6000$ Kuramoto calculations, starting from states with oscillators of fully random phases the growth exponent at the synchronization transition point was found to be $\eta =0.55\left(10\right)$, away from the mean-field expectations by scaling relations: $\eta =0.75$ [17]. It was suspected that finite size scaling and corrections hid the true mean-field behavior numerically in [16] as well as in [17]. Here we provide a more extended numerical study with the aim of clarifying this .

2. Materials and Methods

The model introduced by Kuramoto [5] is one of the most studied for oscillatory systems. Besides the original global coupled system one can define a locally coupled version on graphs [3], in which phases ${\theta }_i\left(t\right)$, located at the N nodes of a network, follow the dynamical equation

$$\dot{{\theta }_i}\left(t\right)={\omega }_i^0+K\sum _j{A}_{ij}sin[{\theta }_j\left(t\right)-{\theta }_i\left(t\right)].$$

(1)

The global coupling K is the control parameter, by which one can tune the system between asynchronous and synchronous states. The summation is performed over the nearest neighboring nodes, with connections described by the adjacency matrix $A_{ij}$ and ${\omega }_i^0$ denotes the quenched self frequency of the i-th oscillator, chosen from a Gaussian distribution, with zero mean and unit variance. When $A_{ij}$ describes a full graph, the dynamical behavior is mean-field type [18]. The critical dynamical behavior has been explored on various random graphs [16,17]. In regular lattices synchronization phase transition can happen only above the lower critical dimension $d_l^{-}=4$ [6]. In lower dimensions, a true, singular phase transition in the $N\to \infty$ limit is not possible, but partial synchronization emerges with a smooth crossover if the oscillators are coupled strongly.

To investigate the relaxation to the steady state, we measured

$$z\left(t_k\right)=r\left(t_k\right)exp\left[i\theta \left(t_k\right)\right]=1/N\sum _{j}exp[i{\theta }_j\left(t_k\right)],$$

(2)

where $0\le r\left(t_k\right)\le 1$ gauges the overall coherence and $\theta \left(t_k\right)$ is the average phase at discrete sampling times $t_k$, which was chosen to follow an exponential growth: $t_k=1+1.{08}^k$ to spare memory space.

The sets of equations (1) were solved numerically for ${10}^3-{10}^4$ independent initial conditions, initialized by different ${\omega }_i^0$-s and random initial phases ${\theta }_i\left(0\right)$ were used. Then sample averages for the phases and the frequencies give rise to the Kuramoto order parameter

$$R\left(t_k\right)=\langle r\left(t_k\right)\rangle ,$$

(3)

which is non-zero above the critical coupling strength $K>K_c$, tends to zero for $K<K_c$ as $R\propto \sqrt{1/N}$, or exhibits a growth at $K_c$ as

$$R(t,N)=N^{-1/2}{t}^{\eta }{f}_{\uparrow }(t/N^{\tilde{z}}),$$

(4)

in case of an incoherent initial state, with the dynamical exponents $\tilde{z}$ and $\eta$. We also measured the variance of the frequencies

$$\mathsf{\Omega }\left(t\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N{(\overline{\omega }\left(t\right)-{\omega }_j\left(t\right))}^2,$$

(5)

which is expected to exhibit an asymptotic decay law $\mathsf{\Omega }\left(t\right)\propto t^{-d/2}$ in the ordered phase, by linear approximation [17], and is also confirmed numerically in the nonlinear regime [17,19].

The Kuramoto differential equation system (1) was solved using the adaptive Bulirsch–Stoer stepper [20] implemented with support for CUDA GPUs boost::odeint [21] via the VexCL library [22]. Node and edge information of sparse random graphs were stored using a memory layout optimized for the efficient parallel evaluation of the Kuramoto computation of interaction term. This implementation allowed the numerical treatment of large systems with $N=12{000}^2$ nodes and containing up to $\sim 6·{10}^8$ edges, while sampling thousands of realizations per graph type on a small GPU cluster.

2.1. The power-Law Decaying Long Range Graph (ll)

Graphs with variable topological and spectral dimension can be generated by the addition of long-range connections with probability decaying with the distance l algebraically as

$$p_l\propto \beta l^{-s},$$

(6)

parametrized by the exponent s, see for example [23]. We used this method by the addition of links to two dimensional lattices, with periodic boundary conditions and linear sizes $L=12000$. For $s<4$, infinite graph dimensions are obtained. The probability of long-range links decreases as s increases and for $s>4$, any long-range links are suppressed to simply give the original underlying two-dimensional lattice. This construction provides the possibility to generate systems with their dimensionalities vary from 2 to ∞ at $s=4$ by tuning $\beta$. As for the spectral dimension the situation is different as we will show it in Sect. 3.1

2.2. Spectral Analysis

The spectrum of the Laplacian matrix of a complex network encapsulate vital information about the network’s structural properties, such as its connectedness and connectivity, its ability to synchronize or conduct diffusion processes, and its resilience to structural perturbations, etc [24]. In particular it is related to the linearized Kuramoto equation, which describes a random walk. Such spectra were also shown to encode the dimensionality information of a network [25]. Following Refs. [11,14], we adopt the normalized Laplacian L with elements

$$L_{ij}={\delta }_{ij}-A_{ij}/k_i$$

(7)

for unweighted networks, where $k_i$ denotes the degree of node i. The normalized Laplacian has real eigenvalues $0={\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_N$ that form a spectrum, the density of which is characterized by the following scaling behavior [11,25]

$$\rho \left(\lambda \right)\simeq {\lambda }^{d_s/2-1}$$

(8)

for $\lambda \ll 1$, where $d_s$ is the spectral dimension. The cumulative density is then given by

$${\rho }_c\left(\lambda \right)={\int }_0^{\lambda }d{\lambda }^{\prime }\rho \left({\lambda }^{\prime }\right)\simeq {\lambda }^{d_s/2}.$$

(9)

Note, that the smallest nonzero eigenvalue ${\lambda }_2$ or the Fiedler value quantifies the connectivity of the network. In connected networks with finite spectral dimensions, the dimension $d_s$ is an overall measure for the local connectivity, while as the system size N increases, the resulting larger network will also have more room for disconnection, giving rise to a smaller ${\lambda }_2$. Therefore, ${\lambda }_2$ depends on both $d_s$ and N. By imposing ${\rho }_c\left({\lambda }_2\right)=1/N$ (the cumulative density of this smallest nonzero eigenvalue over N eigenvalues), ${\lambda }_2$ then scales with the network size according to the following power law [14].

$${\lambda }_2\sim N^{-2/d_s}.$$

(10)

More generally, in order for the spectrum to display the power-law behavior (8) and (9), each eigenvalue should follow the similar leading-order scaling [26]

$${\lambda }_i\sim N^{-2/d_s}.$$

(11)

3. Results

3.1. Spectral Analysis Results

Equations (8)-(11) all suffice for the estimation of the spectral dimension. However, since $d_s$ values obtained from Equations (8) and (9) may depend on how many low-lying eigenvalues are chosen, these two scaling relations are not very reliable. In practice, it is more apt to utilize the finite-size scaling relations (10) and (11) to obtain $d_s$ by observing how the scaling exponents converge as $N\to \infty$. In Figure 1, we applied Equation (10) to estimate the dependence of $d_s$ on N for different s values, which shows that the scaling form Equation (10) only gives a consistent, reliable $d_s$ value by extrapolating to the thermodynamic limit $N\to \infty$. It also shows that the intended upper critical dimension $d_{sc}=4$ for this study corresponds to choosing $s\approx 3$.

To more systematically determine how $d_s$ varies with s in the thermodynamic limit, we may estimate $d_s$ according to Equation (11) with respect to several eigenvalues for large enough system sizes. From the results of Figure 1, taking $N={140}^2,{160}^2,{180}^2$, and ${200}^2$ should give reasonable estimations, while pushing the lateral size L to even larger values will be rather computing resource demanding for eigenvalue calculations. The resulting Figure 2 (see the caption for the calculation details) shows that $d_s$ decreases monotonously with respect to s and tends to the dimensionality of the underlying two-dimensional lattice for large s, as expected according to Section 2.1. This figure further justifies that $d_s\approx 4$ for $s=3$ and $d_s\approx 2$ for $s=4$ in agreement with the relevancy of the long links. In Section 3.3, we shall take $s\le 3$ as we intend to study synchronization dynamics on networks with $d_s\gtrsim 4$.

3.2. Regular Lattice Kuramoto Order Parameter Behavior

First we solved Equation (1) in case of 5 dimensional lattices, where the mean-field type of behavior is expected. We initialized the solver by random initial phases and followed the evolution of $R\left(t\right)$ up to $t_{max}=4000$ time steps. We performed these calculations for $L=32,40$ linear sizes, for thousands of independent random initial conditions. Here we show the results of $L=40$ on Figure 3, because for $L=32$ the cutoff happens at earlier times. Using a local slope analysis, defined by the logarithmic derivatives of the growth Equation (4) at the discrete time-steps $t_k$, near the transition point

$${\eta }_{\mathrm{eff}}={\displaystyle \frac{lnR\left(t_{k+4}\right)-lnR\left(t_k\right)}{ln\left(t_{k+4}\right)-ln\left(t_k\right)}},$$

(12)

we estimated $K_c$ as well as the exponent $\eta$. After selecting a proper rescaling of the horizontal axis, we can observe a linear infelxion curve at $K=K_c=2.046\left(1\right)$, separating the up and down veering curves, corresponding to super and sub-critical cases and we can extrapolate to $\eta =0.77\left(2\right)$ in the $t\to \infty$ limit. The initial growth $\propto t^{1/2}$ behavior, similar as in [17], leads to a strong corrections to scaling. Note, that simple power-law fit on the $R\left(t\right)$ would provide a smaller growth exponent. Our results is the first one to show clearly the expected asymptotic mean-field scaling, because analyses on smaller sizes hindered to see it due to the corrections and early time cutoff.

We repeated this analysis for $d=4$, but here the expected logarithmic corrections at $d=d_c$ make the scaling even more complex. We show only the effective exponents for the calculations of ${100}^4$ lattices on Figure 4, again after a sample average of 1000-5000 independent initial conditions. By rescaling the horizontal axis on the logarithmic scale we can see an inflexion curve at $K\simeq 4$ that may indicate the expected behavior at $d=d_c$, but care must be taken due to unknown finite size cutoff. Thus we just claim that our numerical results canconfirm the mean-field behavior, with logarithmic corrections.

3.3. The ll Graph Order Parameter Results

Now we show our results for the ll graph, where we used large two dimensional base lattices, with maximum linear sizes $L=12000$ and added PL decaying random long links with exponent $s=3$, corresponding to $d_s\ge 4$, according to the spectral analysis shown in Section 3.1. Thus, one would expect to see a simple mean-field behavior, instead as Figure 5 shows, we find a smeared transition, without any clear separatrix extrapolating to the mean-field value $\eta =0.75$ in the infinite time limit. Unlike in the regular lattice cases the local slope curves seem to saturate in the region $0.85<K<0.98$ and $0.3<1/t<0.03$ before the finite time cutoff. But the curves also exhibit wavy modulations, which persist even after averaging over 2000 independent samples, and which prohibits to extend the calculations further.

Such wavy modulations, non-monotonic corrections in the local exponents were also present in case of Kuramoto on random Erdős-Rényi graphs [9], where the authors could not clearly justify relevant deviations from the mean-field behavior, which is present here. Solving Kuramoto on ll graphs with smaller s values runs into rapidly increasing computational difficulties, as it requires larger memory and the communication via long links becomes slow.

We have also obtained results for $s=2.9$ resulting in a higher spectral dimension, where the deviations from the mean-field are even more obvious as shown on Figure 6.

Although we have results for a limited number of couplings, and the corrections are wavy, the results suggest a critical point at $K_c\approx 1.8$, with the corresponding exponent $\eta =1.18\left(3\right)$, far away from the mean-field value ${\eta }_{MF}=0.75$. This resolution of data does not allow us to see a stretched (smeared) dynamical critical region as in case of $s=3$.

For $s=2.8$ the scaling region of $R\left(t\right)$ is even more narrow, instead of that we show in the inset of Figure 6 the evolution of the frequency spread, multiplied by $t^{2.6}$ in system with with $L=6000$, for which PL decay can be found in case of $K=1.6$. The linear approximation predicts $\mathsf{\Omega }\left(t\right)\propto t^{-d/2}$ for Euclidean lattices, below $d_c$ [17], and the observed behavior agrees well with the spectral dimension $d_s\approx 5$, in case of $s=2.8$ and contradicts with the mean-field behavior $\mathsf{\Omega }\left(t\right)\propto t^{-2}$ again.

4. Discussion

In conclusion we have compared the synchronization transition of the Kuramoto model at and above the upper critical dimension and found numerical evidence that the quenched topological disorder is relevant and a smeared, frustrated phase transition seems to emerge instead of the mean-field behavior. This appears even more clearly on lattices of the $d_s\ge 4$ spectral dimensions. Note, however that an oscillating correction to scaling is also present in the effective exponents, which makes it more difficult to analyze data even if we consider very large sizes and averages over thousands of independent initial conditions. Previously such oscillations were also shown in earlier numerical results of the dynamical Kuramoto [17] model on lattices, as well as on random graphs [9]. To show the dimensional equivalence we determined the spectral dimension on ll graphs, but note that topological and graph dimensions are expected to be equal on regular graphs. The graph dimension of the $s\le 3$ ll graphs considered here is infinite. The existence of such frustrated phase transition at high dimension is crucial for understanding models of the brain and warrants further investigations in this direction.