Explosive synchronization coexists with classical synchronization in the Kuramoto model

Synchronization is a ubiquitous feature in a wide range of complex systems. Many physical systems exhibit synchronization: fireflies flashing or people clapping in unison,^1,2 diurnal rhythms in human organ systems, and AC electricity generation,^3–5 to name but a few. Several canonical models have been proposed to capture the way in which complex systems composed of many entities spontaneously synchronize. One of the most widely studied is the Kuramoto oscillator model,⁶ where each oscillator fulfills

in which θ_i is the instantaneous phase, ω_i is the natural frequency of oscillator i, λ is a parameter representing the coupling strength, and A_ij is an adjacency matrix describing the topology of the connections between the nodes. In this way, all nodes simultaneously pull their neighbors toward their phase and are themselves pulled towards their neighbors' phase. For many network topologies, as the coupling strength λ is increased, the oscillators spontaneously synchronize to the same frequency.^7,8 Recently, explosive synchronization has been observed in coupled networks. In contrast to classical synchronization transitions, explosive synchronization (ES) is discontinuous and irreversible. Since abrupt synchronization transitions (albeit without hysteresis) were first reported in 2005,⁹ explosive synchronization has been studied widely.^10–21 For instance, it was studied in the context of periodic phase oscillators for scale-free (SF) network's topologies, where an imposed positive correlation between the natural frequencies of the oscillators and the degrees of nodes was shown to lead to the explosive transition,¹⁰ which was later verified experimentally in a chaotic electronic system.¹¹ A similar transition was observed when natural frequencies and coupling strengths were correlated.^18,19 Both correlation approaches can be shown to be equivalent in the mean-field.^17,19 Furthermore, it was also shown that ES can be considered as the counterpart of an explosive percolation process^22,23 in dynamical phase space.²⁴ 

Another major direction of research on complex systems in recent years is the expansion of network science^25–30 to the study of multilayer networks and networks of networks.^31–35 Essentially, a network of networks is a system in which a set of nodes are connected via qualitatively different links. So far, the combination of links that has been most fruitful to study has been percolation on a network composed of connectivity and dependency links. The connectivity links represent the generic flow that is assumed in single network studies of percolation^29,36 while the dependency links represent the supply of some critical resource, without which the target node is unable to function. It was found that this combination of links can lead to cascading failures and abrupt, first-order transitions, which are absent in single network percolation. This framework was first introduced in 2010 to describe the dependencies between different infrastructure networks, where the connectivity links are exclusively within the network while the dependency links connect between the networks.³¹ However, it was later shown that the same effects arise in systems where the dependency and connectivity links are combined in a single network.^37–39 

One of the most promising models that exhibits explosive synchronization phenomenon is the adaptive coupling introduced by Zhang et al.⁴⁰ In this model, the coupling strength of a given node is adaptive: it changes as a function of the dynamic state of the system. Specifically, the standard coupling λ is multiplied by a factor α_i which is proportionate to the local-synchronization state—either at node i or at some other node. This model is significant because it requires no special constraints on the frequencies or topologies but is in fact totally generic. This adaptive coupling represents an auto-catalytic interaction if it is to the same node: as a node's neighborhood becomes more synchronized, the node itself becomes easier to synchronize. Or if the link is to another node, it represents an interdependent or excitatory relationship (in the case of neural models) between the nodes. It is of particular interest because, by tying the coupling strength of one node to the local-order of another node, it effectively introduces a new type of link, and can be understood in the framework of networks of networks.⁴¹ 

Because the abrupt, first-order, transitions, which characterize explosive synchronization in the adaptive coupling model of Zhang et al.⁴⁰ are absent when the coupling is removed, it is of fundamental interest to understand at what point the transition between the smooth transition behavior and the abrupt transition takes place. The analogous question has been studied extensively for percolation in interdependent networks, where it has been shown analytically that for random networks there exists a critical coupling fraction below which the system undergoes a continuous second-order transition and above which it undergoes an abrupt first-order transition.^42–44 Spatially embedded networks, on the other hand, undergo a first-order transition for any finite level of coupling when the dependency links are random,⁴⁵ but have a critical coupling fraction that depends on the length of the dependency links when those are of limited length.^46,47 Similar studies have been conducted for explosive synchronization, and it has been shown that there exists a critical fraction of nodes with correlated degree and natural frequency in scale-free networks below which the transition is not explosive.¹² 

Zhang et al.⁴⁰ analyzed this model for full coupling (f = 1) but left open the question of whether or not there exists a minimal coupling f_c below which the explosive transition does not occur or whether the explosive transition co-exists with the standard Kuramoto synchronization transition. Here, we develop a mean-field theory for partially coupled oscillator networks and show via analytic calculations and numeric simulations on very large networks that f_c = 0 and the phases coexist. This is the first result that shows explosive synchronization coexisting with the standard phase-transition with no critical coupling threshold. As such, it has broad implications for the realizability of explosive synchronization in real-world systems. Because the phenomenon emerges for any finite fraction of adaptively coupled nodes and without specific constraints on frequencies or degrees, it is more likely to be observed in natural systems.

We describe a system of coupled oscillators which behave according to the following dynamic equation:

where θ_i and ω_i are the phase and natural frequency of oscillator i, respectively, A_ij = A_ji = 1 if oscillators i and j are connected and zero otherwise, λ is the overall coupling strength and α_i is equal to the local order parameter r_i for a fraction f of the nodes and equal to 1 for the remaining 1 − f. The local order is defined as

where k_i is the degree of node i, ψ is the average phase of the neighbors of i, and the sum is taken over all of the immediate neighbors of node i. The new term, α_i, introduces an adaptive factor to the coupling strength which gives rise to the explosive transition. When α_i = r_i, the onset of synchronization is auto-catalytic: as the neighbors of node i become more synchronized, their ability to impact the phase of node i also increases. In general, the local order for the α_i term can be measured at node i itself or at another node, even a node in another network. In this sense, the α_i term can be considered a link (possibly a self-link) of a qualitatively different sort than the typical Kuramoto oscillator coupling.

For the case of full adaptive coupling (f = 1), Zhang et al.⁴⁰ obtained the following self-consistent equation for single networks or pairs of identical networks based on a mean-field approximation

or in the continuum approximation

Over the domain of R, this equation has one or three solutions depending on the value of λ. In the mean field, the result is the same whether the adaptive term is measured at the node itself, at a different node, or even at a node in another network with the same degree distribution.

Though that work included numerical results for systems with adaptive coupling on only a fraction f of the nodes, they did not treat this case theoretically. Furthermore, the numerical tests there appeared to indicate that below a critical fraction f_c, the explosive transition disappears and the system undergoes continuous onset of synchronization. However, the authors were unable to tell if those results were indicative of a critical threshold f_c or merely large fluctuations in the order parameter that overwhelm the hysteresis region. Here, we extend the mean-field theory for a fraction f of adaptively coupled oscillators and 1 − f standard Kuramoto oscillators. We find that, in fact, the explosive transition occurs for all finite fractions of adaptively coupled nodes. This was not observed in the earlier study because, though present, the explosive regime is very small for f and essentially undetectable in small systems. To overcome this problem, we have shifted from CPU to GPU implementations, and have thus been able to simulate systems up to 1000 times larger (∼10⁶) than the previous research (or other standard studies in the field).

Because the theory is the same for single networks with self-coupling and identical pairs of networks, we present the solution for a single network for simplicity. We assume that there are two populations of oscillators: a fraction f which are adaptively coupled (I) and 1 − f non-adaptively coupled (II). For concreteness, we assume that ω_i are drawn uniformly from a uniform distribution between −1 and 1 and that the network has a Poisson degree distribution. This gives us two equations:

These equations can be rewritten using the local order parameter (2) and average phase of neighbors ψ as

In the mean-field approximation we approximate $ri≈R, ψ≈Ψ$ (the average global phase) and obtain N uncoupled equations:

Changing variables to $Δθi=θ−Ψ$ with $Ψ̇=Ω$⁠, the average natural frequency, and utilizing the fact that the natural frequency distribution is centered at zero (i.e., Ω = 0), we obtain

The $Δθi$ values are self-averaged in the calculation of R

In the thermodynamic limit, only the phase-locked oscillators contribute to this sum, and we determine which oscillators are phase-locked by the fixed-points

This gives the conditions of

for the two populations of oscillators. We can now write a combined self-consistent equation for R with both populations,

where we have used the fact that the fraction f of adaptively coupled oscillators is selected independently of the degree k_i and natural frequency ω_i. In the mean field approximation, we sum over all of the oscillators as if they are of each type and combine the sums with the appropriate weights.

To make the calculation tractable, we take the continuum limit to turn the sums into integrals

In the continuum limit, instead of summing the individual ω_i and k_i, we integrate with the appropriate distributions (⁠$g(ω)$ and P(k), respectively). Because we are treating the case of uniformly distributed frequencies and Poisson degree distribution, we have $g(ω)=1/2$ between −1 and 1, 0 elsewhere and $P(k)=k¯ke−k¯/k!$⁠. First thing we note is that the integrated terms of the internal (dk) integral enter via P(k) or in the combinations

So we simplify Eq. (19) as

where

We note immediately that if $R=0$ then $α=0$ which brings the boundaries of integration to 0 and trivially fulfills the condition. Therefore, R = 0 is always a solution, though it is not necessarily always stable. Solving this integral numerically, as described in the  Appendix, we can find solutions for R and obtain analytic predictions for all values of f, see Fig. 1.

For all values of f, there are between one and three solutions for R, depending on λ as shown in Fig. 1. The one solution which is always present for all f and all λ is R = 0 (the desynchronized phase), as noted above. For all f, we can define a point λ_c for which below λ_c R = 0 is the only solution, but above which the stable phase appears. In the special case of f = 0, only this branch rises from the R = 0 branch. At this point, the zero solution becomes unstable and the only stable solution is the synchronized solution. For all other values of f, there is a bifurcation in the branches, as shown in Fig. 1. Even after the bifurcation point, the desynchronized state remains stable to small fluctuations as long as an unstable solution separates it from the synchronized state. Therefore, as λ is increased to λ_c, the system must transition discontinuously to the new branch. However, this transition does not occur as soon as the second branch appears, because the desynchronized phase remains stable to small fluctuations.

When we use a fixed-point analysis of the solutions, we find that letting the rhs of Eq. (22) be F(R)

for all solutions in the top branch and zero branch, whenever a middle branch exists and they are thus stable to infinitesimal fluctuations. The middle-branch, on the other hand, is always unstable.

Zhang et al.⁴⁰ found that the unstable branch continues to exist for all $λ>λc$⁠. When f is decreased, we find that there is a critical value $f*(k¯)$ of coupling below which the unstable branch ends at some value $λ*$⁠. When $λ>λ*$⁠, the zero solution is not stable at all and the system transitions spontaneously, even in the absence of fluctuations.

From our numerical tests, as well as the results in Ref. 40, it is clear that the desynchronized phase becomes unstable at a certain value λ_f which is substantially lower than $λ*$⁠. Determining where exactly that happens is not possible in the mean-field approximation because the transition is driven by fluctuations, which we have neglected. However, as λ increases, the characteristic size of the fluctuations increases and once they are large enough to reach the unstable branch the system passes over to the synchronized state. Therefore, as reported in Ref. 40, if the system begins in the desynchronized state at λ = 0 and λ is adiabatically increased, at a certain value λ_f the fluctuations are large enough to pass the unstable branch and the system transitions discontinuously to the synchronized state. Likewise, if the system begins in the synchronized state and λ is adiabatically decreased, the system will remain synchronized until a value λ_b (with $λb<λf$⁠), at which point it spontaneously desynchronizes. Based on the simulations of networks of size $∼106$ (three orders of magnitude larger than those which are typically studied), we suggest that the determining factor of the phase transition is the distance (in order-parameter space) from the initial branch (synchronized or unsynchronized) to the unstable branch. This is similar to the assessment of Zou et al.⁴⁹ who analyzed scale-free networks with explosive synchronization due to correlations between degree and natural frequency and found that the forward and backward transitions were due to the size of the basin of stability. In Fig. 2(a), we have plotted $ΔR$⁠, the size of the gap between the R = 0 solution and the unstable branch which describes the necessary fluctuation size for the forward transition, and in Fig. 2(b) we show the analogous figure for the backward transition, based on the distance from the synchronized stable branch to the unstable branch. For the forward transition, we show that the same jump size characterizes the transition for all values of f. The backward transition is more difficult to evaluate because the numerical fluctuations in the critical threshold coincide with a quickly varying region of $ΔR$⁠, but the hypothesis that it is the same fluctuation size is consistent with the mean-field theory and our numerical measurements.

The assumption that there is a critical fluctuation size $Δr$⁠, allows us to make a prediction for the dependence of the hysteresis region size on the fraction of adaptively coupled nodes f. In Fig. 3, we show how this prediction compares to our measurements of $d=λf−λc$⁠. We find that the prediction based on this hypothesis is accurate to within the deviations due to the mean-field approximation, which are $∼10−3$ (cf. Fig. 1(b)).

In conclusion, we have found that in fact there is no critical coupling for explosive synchronization when it is introduced via the adaptive coupling term as in Ref. 40. Even before this study, the model of explosive synchronization presented by Zhang et al.⁴⁰ was compelling because it required no particular assumptions on topologies or frequency distributions. This result, namely, that any finite fraction of adaptive coupling induces an explosive transition that is discontinuous and exhibits hysteresis, establishes the model of adaptive coupling in the Kuramoto model as a particularly effective way to describe discontinuous jumps in synchronization in real-world complex systems.

M.D. thanks the Azrieli Foundation for the award of an Azrieli Fellowship Grant. O.M. and S.K. thank Russian Foundation for Basic Research (Project Nos. 15-02-00624 and 15-52-45003). S.H. acknowledges the MULTIPLEX (No. 317532) EU project, the Israel Science Foundation, ONR Global, and DTRA for the financial support. S.H. also thanks the Nvidia Corporation for a small hardware grant, which was used for this research. We also thank S. Guan and Z. Liu for several helpful discussions on previous work with this model. We also thank Ivan Bonamassa for many helpful comments on the manuscript.

There are in fact two well defined regions of integration, for Equation (23) for $k<1/α$ and $k>1/α$⁠, as shown in Fig. 4, the total integral can be divided as

Now, by changing variables to $ω′=ω/αk$ and thus $dω=αkdω′$ we can simplify the bounds and decouple the terms in the square root to obtain for I₁

For I₂ the boundaries of ω are −1,1 and we get a less straightforward equation

In this way, $F(k¯,α)$ can be calculated numerically and then plugged in to Eq. (22) to solve the solutions for all f and λ on any given topology.