EXPLOSIVE SYNCHRONIZATION IN MONO AND MULTILAYER NETWORKS

. Explosive synchronization, an abrupt transition to a collective coherent state, has been the focus of an extensive research since its ﬁrst obser- vation in scale-free networks with degree-frequency correlations. In this work, we report several scenarios where a ﬁrst-order transition to synchronization occurs driven by the presence of a dependence between dynamics and network structure. Therefore, diﬀerent mechanisms are shown to be able to prevent the formation of a giant synchronization cluster for suﬃcient large values of the coupling constant in both mono and multilayer networks. Using the Ku- ramoto model as a reference, we show how for an arbitrary network topology and frequency distribution, a very general weighting procedure acting on the weight of the links delays the synchronization transition forming independent synchronization clusters which suddenly merge above a critical threshold of the coupling constant. A completely diﬀerent scenario in adaptive and multilayer networks is introduced which gives rise to the emergence of an explosive syn- chronization when a feedback between the dynamics and structure is operating by means of dependence links weighted through the order parameter.


1.
Introduction. Synchronization is one of the most fascinating processes in complex networks' dynamics. The spontaneous order of the network's units into a collective dynamics is an emergent phenomenon whose appearance relies on a delicate interplay between the topological properties of the network and the main features of the dynamical system associated to each graph's node [6,2].
Considered as a phase transition, in a large majority of cases synchronization occurs as a second-order transition involving a continuous and reversible change into order of the macroscopic state. However, it has been recently reported that under certain circumstances it is possible the observation of explosive synchronization In both panels, the legends report the color and symbol codes for the different plotted curves. In (c) and (d), the degree k i that each node achieves after the network growth is completed (upper plots) and the average of the natural frequencies ω j of the neighboring nodes (j ∈ N (i), bottom plots) are reported vs. the node's natural frequency ω i , for k = 100 and frequency gaps γ = 0.0 (c), and γ = 0.4 (d). The red solid line in (d) is a sketch of the theoretical prediction f (ω). Panels (e) and (f) show r (color coded according to the color bar) in the parameter space (λ, γ) for (e) k = 20 and (f) k = 60. The horizontal dashed lines mark the separation between the region of the parameter space where a second-order transition occurs (below the line) and that in which the transition is instead of the first order type (above the line). The yellow striped area delimits the hysteresis region.
(ES), an irreversible and discontinuous transition to the coherent state [4]. Originally, ES was described in all-to-all coupled ensembles of Kuramoto oscillators [11] with specific distributions of frequencies [17,15]. Later on, various kinds of node's degree and frequency oscillator correlations were found to be able to induce ES in networks of periodic and chaotic oscillators [9,14].
Disclosing the mechanisms at the root of these unexpectedly abrupt transitions is of fundamental relevance for a deeper understanding of the networks' structure and dynamics. For instance, recent clinical evidence seems to indicate that ES is actually the mechanism driving the transition from normal to pathological brain behavior during epilepsy, one of the world's most prominent brain disorders [10]. ES has also been reported to appear in the transition to synchrony in power networks [18].
Due to the relevance to practical applications, ES has attracted an extensive attention, and stimulated a large amount of research. In this paper, our aim is to provide some insight on the recent findings regarding the microscopic mechanisms of ES in complex networks composed of either a single or several layers [4].
2. Explosive synchronization in monolayer networks. In our research of this phenomenon we use as a benchmark a complex network of dynamical units whose instantaneous phase evolves in time according to the simple and paradigmatic model to describe periodic oscillators introduced by Kuramoto [11]: where θ i is the phase of the i-th oscillator with natural frequency ω i and λ is the coupling constant. The topology of the network is uniquely defined by the adjacency matrix A = (A ij ), since A ij = 1 if node i is connected with node j and A ij = 0 otherwise. The frequencies are chosen according to a given distribution g(ω). The level of synchronization can be monitored using as an order parameter the value of r = 1 N | N j=1 e iθj (t) | T , with ... T denoting a time average over a conveniently large time span T .
Early works in ES focused on the effect of imposed correlations between the node's degrees and the corresponding natural frequency of each oscillator [9,14], i. e. the node frequencies were tailored such that ω i ∼ k i in Eq. (1). However, we intend to show that ES is by no means restricted to such rather limited case, but it constitutes, instead, a generic feature in synchronization of networked oscillators.
The key microscopic feature introduced by the k − ω relationship is not the correlation itself, but the subsequent consequence that the network acquires frequency disassortativity in the presence of heterogeneous degree distributions, which is the case of the real-world networks. This frequency disassortativity imposed by the correlation has the secondary effect of hampering connections between nodes whose frequencies are close. As a consequence, the clustering of possible synchronization seeds is suppressed, and the synchronization state is forcefully frustrated up to a very high coupling.
Based on such hypothesis, it is possible to test the minimal condition for the transition from unsynchronized to synchronized states to be abrupt. The basic idea is to build the network avoiding that oscillators behave as cores of a clustering process. The practical realization of such a condition, which evokes the Achlioptas process leading to explosive percolation [1], is to explicitly constrain the frequency difference between each node i and the whole set N (i) of oscillators belonging to its neighborhood [12]: (2) where a threshold γ for the frequency gap is set, and the network is then grown by means of a conditional configuration model approach: after having initially distributed the oscillator's frequency from a given frequency distribution g(ω), pairs of nodes are randomly selected, and a connection between the nodes is formed only if they verify the condition (2). The process is repeated until the network acquires a predetermined mean degree k . Then, the resulting adjacency matrix is used to simulate the Kuramoto model (1). As we look for a potential occurrence of hysteresis in the transition, each simulation is repeated twice: a first one where the coupling parameter λ is adiabatically increased without resetting the dynamical state, and another one where, starting from a synchronized state, λ is progressively decreased up to the point coherence is lost. Both processes are named in the following as forward and backward transitions, respectively.
In Fig. 1 we report the results obtained by setting g(ω) as a uniform frequency distribution in the interval [0, 1]. Panels (a) and (b) show the phase synchronization index as a function of the coupling strength. In particular, Fig. 1(a) (resp. (b)) illustrates the case of a fixed mean degree k = 40 (resp. a fixed frequency gap γ = 0.4), and reports the results for the forward and backward simulations at different values of γ (resp. k ). A first important result is the evident first-order character exhibited, in all cases, by those transitions for sufficiently high values of γ.
It is very revealing a close inspection of the microscopic organization of the resulting networks produced by applying the frequency gap condition. We found the spontaneous emergence of a k − ω correlation feature associated with the passage from a second-to a first-order like synchronization transition. While such a correlation was imposed ad hoc in Refs. [9,14], here the condition (2) creates for each oscillator i a frequency barrier around ω i , where links are forbidden. The reached degree k i will be then proportional to the total probability for that oscillator to receive connections from other oscillators in the network, and therefore to 1 − Fig. 1(c) refers to the case γ = 0 in which no degreefrequency correlation is present, while in the upper plot of Fig. 1(d), instead, it is reported the ES case occurring at γ = 0.4 and the corresponding normalized function f (ω) = 1 − ω+γ ω−γ g(ω )dω , with g(ω) = 1 for ω ∈ [0, 1], and g(ω) = 0 elsewhere, which gives evidence of the emergence of a clear V-shaped correlation between the ω i and k i . The obtained pattern is clearly nonlinear, depends on the frequency distribution g(ω) and, in general, it is V-shaped. Finally, Figs. 1(e) and (f) show that the rise of a first-order like phase transition is, indeed, a generic feature in the parameter space. Here we report r in the λ − γ space, for k = 20 and k = 60 respectively. The horizontal dashed lines in panels (e) and (f) mark the values of γ c , separating the two regions where a second-order transition (below the line) and a first-order transition (above the line) occurs. The fulfillment of Eq. (2) leads to an explosive transition for a very wide class of distributions of the oscillators' natural frequencies, as shown in Figs. 2(a)-(b) for a Rayleigh distribution. The condition of Eq. (2) can be relaxed in several ways. For instance, a frequency gap in the network growth can be introduced as |ω i − ω j | > γ, where . indicates the average value over the neighborhood N (i). The results for this local mean field gap condition in a homogeneous frequency distribution are shown in Figs. 2(c)-(d).
The previous findings reveal a key factor pointing out to the essence of the ES, but it imposes a specific architecture for the network. However, based on the acquired knowledge about the microscopic nature of the ES, the study can be further extended to the case of a network with an arbitrary frequency distribution and architecture, for which the only action that an external operator can perform is a weighting procedure on the already existing links [13] by modifying Eq. (1) aṡ where Ω ij = A ij |ω i − ω j | α is the weight factor to be applied to the (existing) link between nodes i, j, being A ij the elements of the adjacency matrix that uniquely defines the network, and α a constant parameter which possibly modulates the weights. The strength of the i-th node (the sum of all its links weights) is then s i = j Ω α ij . Numerical simulations of the system (3) for a Erdös-Renyi (ER) [8] network of size N = 500 are reported in Fig. 3(a), for several g(ω) within the range [0, 1]. For the simplest case of a uniform frequency distribution, the un-weighted network displays a smooth, second-order like transition [squares in Fig. 3(a)], whereas the weighting factor has the effect of inducing ES in the same network, with an associated hysteresis between the forward (solid line) and backward (dashed line) transitions. Such dramatic change in the nature of the transition is independent of g(ω), as long as it remains defined in the same frequency range [0, 1]: we obtain identical results for symmetric (homogeneous, Gaussian, bimodal derived from a Gaussian) and for asymmetric Rayleigh distribution: with σ=0.2, and a Gaussian centered at 0 but with just having the positive half (and normalized consequently): with σ = 0.25, showing how general is the behavior for homogeneous networks.
We move now to rigorously predict the onset (and nature) of the explosive transition, by analytically examining the thermodynamic limit in which N oscillators form an all-to-all connected (clique) graphθ i = ω i + λ N N j=1 |ω i − ω j | sin(θ j − θ i ). Using the following definitions [13], the dynamical equations governing the clique can be transformed intoθ i = ω i + λA i sin(φ i − θ i ), whose stationary solution in the co-rotating frame reads as The definition of A ω and φ ω implies that In the case all oscillators are close to synchronization, one can assume that cos θ(x) ≈ r, and therefore G(ω) rs(ω), where s(ω) is the strength of a node with frequency ω. Using this approximation, Eq. (8) takes the form which is a second order differential equation whose integration gives F (ω) depending on the function s(ω). For instance, in the particular case we take a uniform distribution g(ω) in the interval [−a/2, +a/2], we can solve Eq. 10 using the previous definitions of F (ω) and F (ω) to explicitly obtain the strength s(ω) as the second order polynomial and therefore, the integration of Eq. (10) gives Using the consistency equation given in (9) and that F (ω) = 2g(ω) sin θ(ω), it results that where H(x) := 4 4+π x 1+x 2 + arctan(x) . Figure 3(b) shows the emergence of this predicted parabolic relationship between the node strengths s i and the natural frequencies ω i of the oscillators when the synchronization transits from a smooth to an explosive one.
This analysis allows us also to determine the discontinuous nature of the transition by finding the relationship between r and λ. We use that which is an implicit equation in r. When λr ≥ 2+π 4+π ≈ 0.72, sin θ(x) ≤ 1 for all x, which means that all oscillators are frequency locked and, then, On the other hand, if λr ≤ 2+π 4+π , only those oscillators with frequency in the interval [−ω * , ω * ] are locked. Being H(x) continuous and monotone, it can be inverted to obtain: ω * := a 2 H −1 (λr), Hence, if we introduce the following functions being µ = λr. Solving this implicit equation in µ we can obtain a solution for r, which, graphically, for a given value of the coupling constant λ, corresponds to the intersecting points between the function I(µ) and straight lines of slope 1/λ passing through the origin. The function I(µ) presents an inflection point at 2+π 4+π , where the curve changes its concavity (see Fig. (4)). The consequence is that, depending on λ, there are three different kinds of solutions: for λ small, only the trivial solution r = 0 is found since the straight line and I(µ) only intersect at µ = 0; when λ is such that the slope of the straight line is tangent to I(µ) (i.e., when λ = 1.03, the red dashed line in Fig. 4(a)); and when λ increases above this value, the hysteresis appears since r has now three values, with two stable solutions (r = 0 and r ≈ 1) and an unstable solution (see Fig. 4(b)). The solution r ≈ 1 appears therefore abruptly when the slope of the straight line is tangent to I(0) (i.e., when λ = 1.43 (blue dashed line in Fig. 4(a)), which is the point where the stable solution r = 0 collapses with the unstable one, becoming unstable (see Fig. 4(b)).

Microscopic dynamics in explosive synchronization.
More detailed information about the microscopic dynamics underlying ES can be found measuring the local order parameter, r ij [20]: where T is the length of the time window. Its value is r ij =1 for any two phase-locked oscillators, r ij =0 for all pairs of fully uncorrelated oscillators, and intermediate values for any two partially correlated oscillators. One can assume that the matrix R = (r ij ) represents a dynamical functional weighted matrix, whose comparison with the real structural adjacency matrix A = (a ij ) helps to reveal the deep interdependence between structure and dynamics. In Fig. 5 the values of r ij are reported for a fully connected (a)-(d), and ER (e)-(h) networks, for four representative λ values. We observe that only small synchronized clusters of oscillators exist for λ < λ c , while a giant synchronized cluster shows up immediately after λ c , indicating that the small synchronized clusters suddenly merge together right at λ c . It is clear that below λ c , links are always generated among those nodes whose frequencies are relatively close. As a consequence, separated clusters form for λ up to λ c , where instead (and suddenly) all clusters merge together to melt into a giant one.
This information guided us to consider the use of an effective network whose structure explicitly reflects the interplay between the topology and dynamics of the original system and, in an individual basis, it can be useful to easily identify the nodes which are actually acting as synchronization seeds [16]. We define the effective adjacency matrix as where ∆ω ij = |ω i −ω j | is the frequency detuning, and ∆ω max the maximum possible detuning present in the system in order to guarantee C ij ≥ 0. The role of each node in the synchronization process can be extracted from the network defined by C = (C ij ) through the calculation of the standard eigenvector centrality measure of C. Thus we obtain an effective centrality vector Λ C , whose i-th component quantifies the potential of node i to behave as a seed of synchronization.
In Fig. 6 the comparison between Λ C and its topological counterpart Λ A (the eigenvector centrality extracted from the original adjacency matrix A = (A ij )) is reported. For ER networks (Fig. 6(a)), the distribution of the components of the vector Λ C as a function of the corresponding natural frequencies of nodes shows the existence of many seeds of synchronization with natural frequencies close to the center of g(ω). This allows characterizing the connection between the microscale and the macro-scale of the network in a much better way than Λ A , whose components are instead uniformly distributed. For scale-free (SF) [3] networks ( Fig. 6(b)), the synchronization seeds are the hubs, and therefore Λ C and Λ A provide essentially the same information.
3. Explosive synchronization in adaptive and multiplex networks. In the previous sections, our study on ES has been focused on the case of single-layer, or monolayer, networks with non-evolving interactions. Real complex systems, however, have far more complicated forms of interactions, which points to the fact that   Figure 7. Forward (black line with squares) and backward (red line with circles) synchronization transitions for a single network with N=1000 and f=1. The insets report the dependence on f of the average transition points λ F and λ B for ten realizations.
the hypothesis of a single-layer, static network may actually result in an overestimation (or an underestimation) of the problem under study. A more accurate analysis resorts on multilayer networks, and dynamical processes on top of them are attracting more and more attention (see, for instance, Ref. [5] for a recent and rather comprehensive review).
In this section we show how ES can emerge from an adaptive process without introducing any correlation into the system a priori [19]. Let us begin by considering the following modification to Eq. (1) with an explicit fraction f of the nodes adaptively controlled by a local order parameter. The evolution of each oscillator is then ruled by:θ where the new parameter α i accounts for an adaptive control of nodes. For a fraction f of the nodes randomly chosen, α i = r i is assumed for each of the selected nodes, where r i is the instantaneous local order parameter of the i-th oscillator, defined as: and φ denotes the phase averaged over the ensemble of neighbors. For the remaining fraction 1 − f of the nodes, α i = 1. Obviously, f = 0 returns the traditional Kuramoto model, while f > 0 indicates that a fraction of nodes are adaptively controlled by the local order parameter. Besides, another significant point is the choice of the natural frequencies ω i of oscillators, which are taken from a random homogeneous distribution g(ω) in the range [−1, 1]. The presence of an abrupt transition with an associated hysteretic loop in r is evident in Fig. 7, indicating the occurrence of ES above a critical value f c . The inset of Fig. 7 reports the dependence on f of the corresponding forward and backward transition points λ F and λ B . The next step of our study is, then, moving to some theoretical analysis, in order to grasp the essential ingredients at the origin of the observed scenarios. To that purpose, we consider the case of Fig. 7, with f = 1 as an example. Substituting the definition of r i into Eq. (19) one haṡ whereΨ = Ω is the group angular velocity. In the mean field approach, r i = r. Letting ∆θ i = θ i − Ψ, we get∆θ i = ω i − Ω − λr 2 k i sin(∆θ i ). If |ω i − Ω| < λr 2 k i , the oscillator i becomes locked to the mean field. For symmetric g(ω), then Ω = 0 and therefore the phase-locked oscillators fulfill: from where it is possible to recover the contribution of the locked oscillators to the order parameter where h(k, ω) = P (k)g(ω) is the joint distribution with P (k) being the degree distribution of the network. The solution of Eq. (23) is shown in Fig. 8, where the presence of an unstable middle branch is responsible for the hysteretic loop associated to ES observed in Fig. 7. This analytic result allows us for a better understanding of the ultimate causes of ES, and in particular of the microscopic mechanisms that are at the basis of the arousal of explosiveness in the transition. Notice that if one considers the Kuramoto model for the common second-order phase transition Eq. (1) and develops the same mean-field treatment, one obtains that the formula for the order parameter is identical to Eq. (23), but the superior integration limit is replaced by |ω| ≤ λrk, which results in the following consequence: for the backward evolution in Fig.7, when the system begins in a synchronized state, one has r r 2 ∼ 1. However, for the forward transition, r 0, and thus r 2 r. Therefore, the difference in the integration domains strongly impacts the fraction of oscillators belonging to the main synchronization cluster. In the usual case of a second-order transition, the oscillators with closer natural frequencies will first form small synchronized clusters and then these clusters will gradually grow up and merge with the increase of the coupling strength, until eventually forming a giant cluster. On the contrary, in the present case, the factor r 2 in the integration domain has the effect of actually suppressing the merging of small synchronized clusters. Thus, with the increase of λ, more and more free oscillators will be attracted to each of the distinct clusters, but these clusters are prevented from merging with each other. Eventually, when no more free oscillators are left, a discontinuous and abrupt behavior of r will show up as a consequence of the sudden collapse of all clusters.
Finally, we extend our study to the case of inter-dependent multiplex networks [19], where each node has a one-to-one partner with the same index i whose governing equations areθ where the superscripts 1 and 2 denote the network layers 1 and 2, and α 1 i and α 2 i represent the coupling of the two layers via dependency links. In details, if node i belongs to the fraction f , α 1 i and α 2 i refer to the local order parameters r 2 i and r 1 i , i.e. α 1 i = r 2 i and α 2 i = r 1 i (α 1 i = α 2 i = 1 otherwise). The frequency distribution g(ω 1 i ) is still random uniform and in the range [−1, 1], while g(ω 2 i ) can be either random uniform or Lorentzian distribution. Figure 9 shows the dependency of synchronization levels r 1 and r 2 as a function of the coupling λ. In the simulations, layer 1 is always a random ER network with k 1 = 12, and we choose different topologies and frequency distributions g(ω 2 i ) of the layer 2 (see caption for details). One clearly sees that in all cases ES occurs simultaneously in both layers. The subplots show the dependence of the hysteresis loop width on f in each case. The apparition of ES seems very robust against the difference in topology and frequency distribution, which again indicates that the correlation between k and ω is not the essential condition for the emergence of ES. If f is below the critical value f c , synchronization returns to a second-order transition. More importantly, these observations can be quantitatively verified via the mean-field theory. In the traditional second-order transition, oscillators with close frequency firstly collapse into small synchronization clusters, which gradually converge towards a giant cluster at the critical coupling strength.  Figure 9. Synchronization parameters r 1 and r 2 vs. λ for two inter-dependent networks with N = 10 3 and f = 1. Squares and circles (triangles and stars) denote the forward (backward) transitions, and the insets report the average width d of the hysteresis loop as a function of f . Here layer 1 is a ER network with k = 12. From (a) to (b), layer 2 is ER network with k = 12 and k = 6, and g(w 2 i ) is a random homogeneous distribution in the range [−1, 1]. From (c) to (d), layer 2 is ER network and SF network with k = 12, and g(w 2 i ) is Lorentzian distribution and random homogeneous fashion.
Recently, furthermore, this kind of ES has been shown to coexist with the standard phase of the Kuramoto model in the thermodynamic limit [7].