Inter-layer adaptation induced explosive synchronization in multiplex networks

It is known that intra-layer adaptive coupling among connected oscillators instigates explosive synchronization (ES) in multilayer networks. Taking an altogether different cue in the present work, we consider inter-layer adaptive coupling in a multiplex network of phase oscillators and show that the scheme gives rise to ES with an associated hysteresis irrespective of the network architecture of individual layers. The hysteresis is shaped by the inter-layer coupling strength and the frequency mismatch between the mirror nodes. We provide rigorous mean-field analytical treatment for the measure of global coherence and manifest they are in a good match with respective numerical assessments. Moreover, the analytical predictions provide a complete insight into how adaptive multiplexing suppresses the formation of a giant cluster, eventually giving birth to ES. The study will help in spotlighting the role of multiplexing in the emergence of ES in real-world systems represented by multilayer architecture. Particularly, it is relevant to those systems which have limitations towards change in intra-layer coupling strength.


I. INTRODUCTION
Recently, an irreversible synchronization process, called explosive synchronization [1, 2], in which a group of incoherent dynamical units is abruptly set in collective coherent motion, has drawn much attention of the researchers [3]. The abnormal hypersensitivity of ES can be perilous in many physical and biological circumstances such as abrupt cascading failure of the power-grid [4], breakdown of the internet due to intermittent congestion [5], abrupt attack of epileptic seizures in the human brain [6] and abrupt episodes of chronic pain in the Fibromyalgia human brain [7], to name a few. An anesthetic-induced transition to unconsciousness [8,9] and bistability in Cdc2cyclin B in embryonic cell cycle [10] are a few other instances of an abrupt transition. The occurrence of ES transition has also been demonstrated experimentally [11][12][13]. The emergence of ES is shown to be rooted in inertia [14] and a microscopic correlation between frequency and degree or coupling strength of the networked phase oscillators [1,15]. Recently, Zhang et al. [16] showed that a fraction of adaptively coupled phase oscillators gives birth to ES. Adaptation is an inherent feature in the construction of many complex systems, for instance, adaptation in neuronal synchronization of brain regions apropos learning or memory process [17,18]. In many complex systems, the same set of nodes may have different types of connections among them, having each type influencing the functionality of other types. Hence, an isolated network is an unfit candidate to model such systems. Such systems can be precisely framed by a * Corresponding Author:sarikajalan9@gmail.com multiplex network, where different layers denoting different dynamical processes are interconnected by the same set of nodes denoting interacting dynamical units [19][20][21][22][23][24][25][26][27]. For instance, social networks, neuronal networks, and transport systems in which individuals, neurons, and cities have different types of connections among them forming different layers [28]. The multiplex framework has been remarkably successful in providing insights into the dynamical behavior of various processes such as percolation [29], epidemic [30,31], voter model [32], etc., taking place in a combination within a group, commu-nity or population. Recently, the studies apropos the occurrence of ES have been extended to multilayer networks by employing different methodologies, for instance, adaptive coupling [16], intertwined coupling [33], inertia [34], delayed coupling [35], inhibitory coupling [36], frequency mismatch [37]. A few recent studies on multilayer networks have adopted the adaptive coupling dynamics proposed by Zhang et al. [38][39][40] at the heart of their models. It is reported that the occurrence of ES in a multilayer network is exhibited by a fraction of nodes adaptively coupled via local order parameter, within the multiplexed layers, having virtual inter-layer links [16]. In the current work, taking an altogether different cue, we propose an approach in which inter-layer links of a multiplex network are adaptively coupled through global order parameters of the interacting layers, which brings about ES with an associated hysteresis. Our model is sensitive to frequency mismatch between the interconnected nodes which, along with inter-layer coupling strength, determines the size of emergent hysteresis. Furthermore, the proposed model is robust against a variety of network topology in giving birth to ES. We corroborate our findings by providing rigorous mean-field analysis and a good match between the analytical prediction of the order parameter and its numerical evaluation.

II. MODEL AND TECHNIQUE
We consider a multiplex network consisting of two layers, each one having N nodes represented by Kuramoto oscillators. Each node in a layer is adaptively linked with its counterpart in another layer through multiplexing strength which is a function of a measure of global coherence of the interacting layers. The time-evolution of phases of the nodes in the multiplexed layers a and b is ruled by [41]: (1-a) and ω a(b) i represent phase and natural frequency of i th node in layer a(b), respectively. ij =0 otherwise. λ represents the inter-layer coupling strength or multiplexing strength. The multiplexing strength between the two layers is adaptively controlled by the prod-uct of global order-parameters (r a r b ), a measure of degree of phase coherence among the nodes given by where 0 ≤ r a(b) ≤ 1. r a(b) = 0 corresponds to a random distribution of the nodes over a unit circle, whereas r a(b) = 1 corresponds to exact phase synchronization. Eq.1 is numerically solved using RK4 method with timestep dt = 0.01. To determine phase coherence among nodes of layer a, time average of r a is taken for 10 5 steps after neglecting initial 10 5 steps. Initial phases and natural frequencies of the nodes in layers a and b are drawn from a uniform random distribution in the range ≤ γ, respectively. In this work we take γ = 0.5.

III. RESULTS
To investigate the effect of adaptive coupling between all the pairs of mirror nodes, we consider a multiplex network of two globally coupled networks, otherwise mentioned elsewhere. Furthermore, oscillators in the two layers have identical natural frequency distribution but in general, ω a i = ω b i . Fig.1 depicts behavior of r a − σ profile for the considered multiplex network. In the absence of adaptive multiplexing, the usual static λ gives rise to a continuous phase transition in layer a, ( Fig.1(a)). However, it unfolds that the presence of adaptive multiplexing strikingly leads to a discontinuous phase transition (ES) accompanied by a hysteresis, ( Fig.1(b)). Factors determining hysteresis width: Fig.2, further elaborates on how inter-layer coupling strength (λ) and frequency mismatch between the mirror nodes affects the transition to synchronization in layer a. It turns out that an increase in λ widens the hysteresis width ( Fig.2(a-c)) associated with the emergent ES. Similarly, for a given λ, an increase in frequency mismatch, ∆ω, between mirror nodes also widens the hysteresis width ( Fig.2(d-e)). Only the backward critical coupling (σ b c ) manifests a decrease with an increase in λ as well as ∆ω while the forward critical coupling (σ f c ) remains almost the same. To measure the strength of frequency mismatch between the mirror nodes, we con- To obtain a desired value of ∆ω, starting with ∆ω = 0, two pairs of mirror nodes are chosen randomly and their natural frequencies are swapped within the layers. After each swapping, ∆ω is calculated and the change is accepted if the newer value of ∆ω is closer to the desired value, otherwise, the change is discarded. This process is repeated until we get a desired value of ∆ω.  Fig. 3 it is around 0.2. Robustness of ES against network topology: The employed technique is robust against a variety of topologies such as regular ring network, Erdös-Rényi (ER) random network [42] and scale-free (SF) network [43] chosen for the individual layers. For numerical simulations, we replace intra-layer coupling σ a(b) /N in Eq. 1 by σ a(b) [44]. Fig. 4 (a-c) shows that a multiplex network comprising either ER-ER network, ER-regular ring network or ER-SF network does exhibit ES transition with hysteresis as depicted by globally connected layers in Fig. 1 and Fig. 2. However, hysteresis size may depend on the network topology as found in the case of ER-SF networks in Fig. 4

(c).
Mean-field analysis when ∆ω = 1: To understand the mechanism behind the origin of ES, we analytically derive r a(b) for the synchronized state. The natural frequency of each node i in layer a is drawn from a uniform random distribution from interval ω a i = −0.5 + (i − 1)/(N − 1), i = 1, 2, . . . , N . Now, we take ∆ω = 1, i.e., ω b i = −ω a i , Eq.1 can then be rewritten in terms of order parameter (Eq. 2) as: where i = 1, . . . , N , and A = r a r b . In the synchronous stateθ can observe that r a = r b = r. On assuming ψ a = ψ b = 0 [45], Eq. 3(a) results into After substituting sin(2θ a i ) = ±2sin(θ a i ) 1 − sin 2 (θ a i ), Eq.
4 results in a fourth order polynomial The roots of the polynomial are [46] x where Here z i is a real root of the polynomial z 3 [46]. Out of the four roots, a root with physically accepted solution is selected, as follows: Eq. 4 suggests that for any two nodes i, j if ω a j = −ω a i , we get θ a j = −θ a i . Therefore, for the synchronous solution, we consider only those roots of the polynomial which satisfy this condition for phases. Eq. 6 implies that if ω a j = −ω a i , then R i = R j as q i is independent of ω a i and the third order polynomial also does not depend on the sign of ω a i , hence, z i is also independent of sign of ω a i . Moreover, for any ω a i = 0, R i = 0 as for R i = 0 or z i = q i , the third order polynomial would imply r ′ i = 0, which is not possible unless λ is infinite. Eq. 6 also implies that D i = E j , therefore x i1 = −x i3 and x i2 = −x i4 . Thus, a possible pair of the roots are either x i1 , x i3 or x i2 , x i4 . We find that r values corresponding to x i1 and x i3 do not match with the numerical simulations (see star with solid line in Fig. 5 (b)), hence the synchronous state corresponds to x i2 and x i4 . It is quite apparent that for ω a i = 0, x i2 = x i4 = 0. After plugging the values of phases from Eq. 5 into Eq. 2, Eq. 2 can be rewritten as sum of the contributions from the locked and the drifting oscillators as: For a given σ and λ, natural frequency of the locked oscillator must satisfy the following relation If f (θ a ) be the RHS of Eq. 8, then f (θ a ) has extrema at cos (θ a i ) * ± = −σr 8λA ± σ 2 r 2 64λ 2 A 2 + 0.5 obtained from the roots of df /dθ a = 0. Substituting the condition for extrema in Eq. 8, we get |ω a i | ≤ 1 − cos 2 (θ a i ) * ± 3σr ± √ σ 2 r 2 + 32λ 2 A 2 4 (9) RHS in Eq. 9 is larger if we consider cos (θ a i ) * + in the RHS We, therefore, proceed with this choice to consider all possible oscillators into the locked state. In the infinite size limit, the contribution of the drifting oscillators to r is 0, which implies that the second summation in the RHS can be neglected for a large network, and therefore r a ≃ r a locked . In the infinite size limit, the probability of finding an oscillator in layer a at phase θ a , while its mirror node's phase is θ b , is C(ω a )/|θ a |, where is C(ω a ) is a constant [47]. Thus, the second term in RHS can be written as Since integration over ω a can be split up into two parts: −0.5 to −γ c and γ c to 0.5, where γ c is equal to RHS of Eq. 9. For a symmetric natural frequency distribution, g(−ω a ) = g(ω a ), we get r a drif t = 0 from the same arguments as shown in [47]. For a given σ and λ, neglecting the contribution of the drifting oscillators to r a , we find the roots (r a values) of Eq. 7 numerically. Note that for a symmetric natural frequency distribution, one can use either x i2 or x i4 for ±ω a i in Eq. 7. Fig. 5(a, b) demonstrates that the analytical prediction of r a is in a fair agreement with its numerical estimation. In Fig. 5(b), solid and dashed lines correspond to x i 2(4) , while the line with stars corresponds to x i 1(3) . If there exist two r values for a given σ and λ, the line with stars shows only the largest of them. It is obvious that only x i 2(4) represents the synchronous state. Fig. 5(a) depicts that in the absence of adaptive multiplexing, r a gradually increases yielding a partially synchronized state, a single cluster of all the nodes havingθ a i = 0, which grows in size with recruiting more and more nodes in it as σ increases. In the presence of adaptive coupling, Eq. 7 does not have any non-zero solution for σ < σ b c , whereas at σ b c one observes an abrupt transition to r a ≈ 1 (see Fig. 5(b)). In Fig. 5(a), the difference between the numerics and the analytical solution corresponding to σ ≈ 0.4 is approximately 0.1, which might be arising due to the omission of the drifting oscillators in Eq. 7. However, for smaller and larger values of σ, the difference is negligible. The solid line in Fig. 5(b) refers to a stable synchronized state, and the dashed line joining the stable state to the incoherent state, refers to an unstable state [15,16,48]. For any given value of σ, r a = r = 0 is also a solution of Eq. 7 although not shown in Fig. 5(b). At r = a = 0, RHS of Eq. 9 is zero, hence r = 0 is a solution of Eq. 7. The presence of the unstable state is an indicator of simultaneous presence of two local attractors, i.e., incoherent state and the coherent state [16,49]. The values obtained for σ f c from numerical simulations presented in Fig. 1 (b), Fig. 2 and Fig. 5(b) can also be perceived from [50], which shows that the incoherent state loses its stability at σ f c = 4γ/π = 0.636 in the thermodynamic limit. We find that adaptive coupling suppresses the formation of the giant cluster, which in turn, leads to ES. To demonstrate the suppression, we evaluate RHS in Eq. 9 for the two cases, A = 1 and A = r 2 . For the forward continuation of σ, we start with r ≈ 0, and therefore we can safely consider r ≈ 0. Further to simplify the calculations we take σ = λ. For A = 1 and r 2 , cos (θ a i ) * + ≈ 0.7 and 0, and hence sin (θ a i ) are almost the same in both cases. However, the bracketed term in the RHS for the two cases is approximately √ 2λ, 0, and therefore, the adaptive coupling suppresses the number of nodes in the partial synchronous state or in the other words, formation of the giant cluster is suppressed until σ f c is reached. Furthermore, Eq. 9 exhibits that an increase in λ leads to an increase in the RHS of the Eq. 9, indicating that the same number of nodes can get synchronized even at a smaller σ value, which in turn causes a decrease in σ b c . For both the layers being identical, a decrease in the hysteresis width with a decrease in ∆ω can be understood from the synchronous state. Substitutingθ a i =θ b i = 0 in Eq. 1 yields ∆ω = 0, θ a i = θ b i , leading to the cancellation of the inter-layer coupling terms and the two layers become isolated networks, which, in the thermodynamic limit manifests a discontinuous phase transition without hysteresis [51].
Mean field analysis when ∆ω = 1: Using the mean field analysis, we prove that there exists a decrease in σ b c with an increase in ∆ω, as well as there exists a discontinuity in r a for globally connected layers. Except at the higher λ values, in the synchronized state r a 0.8 ( Fig. 2(a-c)), indicating existence of either a global synchronized state or almost absence of drifting oscillators. The phases in this state can be written from Eq. 14 (AP-PENDIX A). We solve Eq. 7 numerically by substituting sin(ψ a − θ a i ) values from Eq. 14. While solving Eq. 7, only real terms, arising due to the locked oscillators, are considered in the summation. Fig. 6 reflects that the numerical and the analytic prediction match very well. Figs. 6(b),(c) show that σ b c decreases with an increase in ∆ω. Furthermore, Eq. 7 does not has any non-zero solu-tion below σ = σ b c ≈ 0.45, (Fig. 6(c)), indicating absence of oscillators havingθ a i =θ b i = 0. The minimum value of r a we get is approximately 0.8, indicating suppression of the giant cluster due to the adaptive coupling. For σ > σ b c , from the two solutions of Eq. 7, only the largest r a value is plotted in Fig. 6 (b),(c) i.e. we ignore the unstable solution depicted in Fig. 5 (b). Note that in the synchronous state, we considerθ a i =θ b i = 0. However, it may be possible that only one of them is zero and other one takes a non-zero value (θ a i = 0 andθ b i = 0 or vice-versa), therefore, contributions from these oscillators have been neglected in the summation in Eq.7. For σ < σ b c , r a ≈ 0 in the numerical simulations suggests negligible number of locked oscillators, while for σ > σ b c , almost all the oscillators are in the locked state. Therefore avoiding the possibility ofθ a i = 0 andθ b i = 0 is a valid assumption. Robustness of ES against change in adaptive scheme: Finally, we demonstrate that ES exists even if we replace the inter-layer multiplication factor r a r b in Eq.(1-a) and Eq.(1-b) by r a and r b , respectively ( Fig. 6(d)). In the infinite size limit, the addition of coupling term proportional to r 2 does not yield any change in the critical coupling at which incoherent state becomes unstable [50]. Therefore, r a r b term in the inter-layer coupling only helps us in fixing σ f c at a constant value, which, otherwise might be sensitive to the parameters λ or ∆ω. Note that in the Zhang et al. model, intra-layer coupling term contains r 2 , which was shown to cause ES in adaptively coupled networks [16], a condition which is not required for the case of adaptive inter-layer coupling.

CONCLUSION
It is known that in a system of networked oscillators whatever microscopic strategy can suppress the synchronization, gives rise to ES. It was earlier reported that a fraction of adaptively intra-linked nodes through local order parameters brought about ES in a multilayer network with dependency links [16]. In the current work, an adaptive inter-linked set up between the pairs of the mirror (interconnected) nodes by means of global order parameter trigger ES in the multiplexed layers. We have discussed in detail how the inter-layer coupling strength, as well as the inter-layer pairwise frequency correlation, enable us to shape the hysteresis of the emergent ES. We provide mean-field analytical treatment to ground the perceived outcomes and substantiate that the analytical predictions are in a fair agreement with the numerical estimations. Our model has some level of similarity with several realworld complex systems having multilayer underlying network structure. For example, large-scale brain multilayer networks can be defined based on the functional interdependence of the brain regions [52]. We hope that our investigation of ES originating from the inter-layer adaptive coupling between layers of a multiplex network would help in advance the understanding of microscopic as well as macroscopic dynamics of adaptive mechanism taking place between intertwined real-world dynamical processes.

APPENDIX A
In general case of ∆ω = 1, one obtains the synchronized layers in which pairs of the mirror nodes are mutually synchronized. For globally connected layers, rewriting Eq. 1 in terms of order parameter (Eq.2): where i = 1, . . . , N . In the synchronous state, we havė θ a(b) i =ω [16], whereω is the mean of the natural frequencies of the entire multiplex network. Furthermore, consideringω = 0 yields synchronous states with a fixed point solution. Numerical simulations suggest that one can assume ψ a ≈ ψ b = ψ and r a ≈ r b = r in the synchronous state ( Fig. 6 (a-c)). While summing up Eq.(11a) and Eq.(11b) taking into account these approximations, yields