A Universal Route to Explosive Phenomena
aa r X i v : . [ n li n . AO ] F e b A Universal Route to Explosive Phenomena
Christian Kuehn a,b and Christian Bick c,d,e a Faculty of Mathematics, Technical University of Munich, Boltzmannstr. 3, 85748 Garching b. M¨unchen, Germany b Complexity Science Hub Vienna, Josefst¨adter Str. 39, 1080 Vienna, Austria c Centre for Systems Dynamics and Control and Department of Mathematics, University of Exeter, Exeter EX4 4QF, UK d Institute for Advanced Study, Technical University of Munich, Lichtenbergstr. 2, 85748 Garching b. M¨unchen, Germany e Oxford Centre for Industrial and Applied Mathematics, Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK (Dated: February 26, 2020)Phase transitions are observed in many physical systems. This includes the onset synchronization in a networkof coupled oscillators or the emergence an epidemic state within a population. “Explosive” first-order transitionshave caught particular attention in a variety of systems when classical models are generalized by incorporatingadditional effects. Here we give a mathematical argument that the emergence of such explosive phenomena isnot surprising but rather a universally expected effect: Varying a classical model along a generic two-parameterfamily must lead to a change of the criticality. To illustrate our framework, we give three explicit examples ofthe effect: in a model of adaptive epidemic dynamics, for a generalization of the Kuramoto model, and for apercolation transition.
Many nonlinear physical systems—ranging from epidemicspreading, synchronization of coupled oscillators to perco-lation on a network—undergo phase transitions as a systemparameter is varied. These transitions can be continuous(second-order) at the transition point or discontinuous (first-order). Discontinuous first-order transitions typically lead toan “explosive” change of system properties; see [1] and ref-erences therein. For a wide variety of systems it has beenobserved that a variation of the model via additional fea-tures leads to the change from a continuous second-order toa discontinuous first-order phase transition. The emergingparadigm is as follows. First, many studies add an addi-tional effect to a classical model. Second, these studies ob-serve that a fundamental change in the phase transition (orbifurcation) structure occurs: Upon variation of a new pa-rameter a previously second-order/soft transition becomes afirst-order/hard transition. As an example consider the classi-cal Kuramoto model, which shows a continuous synchroniza-tion transition. However, varying the distribution of intrinsicfrequencies [2, 3] or generalizing the network to simplicialor higher-order coupling [4, 5] allows for discontinuous syn-chronization transitions. Similarly, adding adaptation [6] orhigher-order coupling structures [7] to models of epidemicspreading can induce a discontinuous transition to the epi-demic state.In this paper, we give a mathematical argument that a tran-sition from a continuous to a discontinuous phase transition isnot surprising but a generically/universally expected effect ifadditional parameters are varied. Specifically, we show thata typical model variation along a two-parameter family of aclassical model with a second-order transition must lead to achange of the criticality to first-order. To illustrate our results,we then demonstrate this effect in three explicit examples ofphysical systems: adaptive epidemic dynamics, synchroniza-tion in the Kuramoto model with non-additive higher-orderinteractions, and a model from percolation theory.
A universal mechanism that modulates transitions.
Wefocus here on dynamical systems that have a mean-field or continuum limit model described near the phase transition byan ordinary differential equation (ODE) x ′ := d x d t = F ( x, y ) , x (0) = x , (1)where x = x ( t ) ∈ R n is the unknown and y ∈ R m are param-eters. Suppose x ∗ = x ∗ ( y ) is a smooth family of equilibriumpoints parametrized by y . For all models we have in mind, onetrivial branch of solutions exists for all parameters so we mayassume upon translation that x ∗ = 0 := (0 , , . . . , ⊤ ∈ R n ,i.e., F (0 , y ) = 0 for all y ∈ R m . Furthermore, supposethat we have a bifurcation point [8] upon parameter variationgenerically given by a single eigenvalue of the Jacobian A ( y ) := D x F (0 , y ) ∈ R d × d crossing the imaginary axis. Using a translation in parameterspace and center manifolds [8] or Lyapunov–Schmidt reduc-tion [9], we may assume without loss of generality that themain bifurcation parameter is p = y and the one-dimensionalODE on the center manifold is given by x ′ = f ( x, p ) , x ∈ R , p ∈ R , (2)with bifurcation point at p = 0 . Then it is well-known that thetwo typical bifurcation points encountered in applications arethe transcritical bifurcation with local normal form x ′ = px + ax , (3)where a = ± determines whether the bifurcation/transitionupon varying p is second-order ( x ≥ , a = − ) or first-order ( x ≥ , a = +1 ). Similarly, if there is an equiv-ariance given by a Z -reflection symmetry in the model via f ( x, p ) = − f ( − x, p ) , then the generic transition is a pitch-fork bifurcation x ′ = px + ax . (4)The pitchfork is second-order if it is supercritical and a = − ,while it is first-order if it is subcritical and a = +1 . pqx FIG. 1. Sketch for the variation of a transcritical bifurcation for thephase space { x ≥ } and parameters ( p, q ) with primary parameter p and second generic unfolding parameter q . Dashed lines indicateinstability of the equilibrium and solid lines indicate stability. Thegrey cases are first-order (explosive, subcritical) transitions, whilethe black diagrams are second-order (non-explosive, supercritical)transitions. A generic model variation with at least one additional freeparameter leads to a vector field f that allows for a change incriticality. Take a generic single-eigenvalue crossing and withthe phase space { x ≥ } . Upon variation of the model, thepersistence of a single-eigenvalue crossing is generic withinone-parameter families of the vector field f . Hence, withoutloss of generality suppose that the single eigenvalue crosses at p = 0 . Furthermore, if we vary the model at least one addi-tional free parameter, say q = y , generically appears. Thisparameter takes into account the additional effect for eachmodel as indicated above. A Taylor expansion at the bifur-cation point now yields f ( x, p, q ) = M X j,k,l =0 c jkl x j p k q l + O ( M + 1) , where O ( M + 1) denotes terms of order M + 1 . The coef-ficients c jkl are constrained: The existence of a trivial branchof equilibria, f (0 , p, q ) = 0 , implies c kl = 0 for all j, k ∈ N = N ∪ { } . Since a single eigenvalue crosses at p = 0 ,we must have ∂ x f (0 , , q ) = 0 , where ∂ x denotes the partialderivative. Hence, we have c l = 0 for all l ∈ N and thus f ( x, p, q ) = c xp + c x + X j + k + l =3 c jkl x j p k q l + O (4) . Now we have all the bifurcation conditions taken care of,one may use bifurcation theory to unfold the singular pointinto a generic family. In particular, the next derivatives of thevector field at the bifurcation point should not vanish. Hence,we must have c = 0 from above and for all other combi-nations of indices that c jkl = 0 if j ≥ and j + k + l = 3 ,where the leading-order non-vanishing conditions are ∂ xxp f (0) = 0 , ∂ xxq f (0) = 0 . This yields the truncated lowest-order two-parameter unfold-ing normal form f ( x, p, q ) = c xp + ( c + c p + c q ) x . (5) pqx FIG. 2. Sketch for the variation of a pitchfork bifurcation for thephase space { x ≥ } and parameters ( p, q ) with primary parameter p and second generic unfolding parameter q . Dashed lines indicateinstability of the equilibrium and solid lines indicate stability. Thegrey cases are first-order (explosive, subcritical) transitions, whilethe black diagrams are second-order (non-explosive, supercritical)transitions. We now apply a scaling (or geometric desingularization, orrenormalization) with a small parameter ε > through thetransformation x xε α , p pε β , q qε γ . For the transcritical normal form (5) we choose α = 1 , β = − , γ = − to obtain (upon a suitable time rescaling) f ( x, p, q ) = c xp + ( c ε + c pε + c q ) x . (6)Hence, one easily checks that there is a sign change of ∂ xx f (0 , p, q ) upon varying q in an interval [ − q , q ] for some q > as long as c = 0 , which we expect generi-cally. Even if c = 0 , we can expand to higher order in q and may thereby eventually change the sign of ∂ xx f (0 , p, q ) so only certain situations with e.g. symmetries and/or non-generic smooth functions could lead to the preservation of thesign for all q ∈ R . Once the sign of ∂ xx f (0 , p, q ) changes,this implies that generically the second parameter is able tochange the transition from second to first order or vice versa.Of course, from the viewpoint of the geometry of the bifurca-tion diagram, this is quite intuitive as shown in Fig. 1 that asecond generic parameter may change criticality.The situation for the pitchfork works very similarly exceptthat an additional symmetry f ( x, p, q ) = − f ( − x, p, q ) has tobe respected. This further constrains the coefficients of theTaylor expansion. Note that if this symmetry is broken thenwe are in the transcritical case if there is still a trivial branchfor all values of the parameters. Hence, we now assume thatthe symmetry holds. Taylor expansion as above gives for abifurcation point with a single eigenvalue crossing f ( x, p, q ) = c xp + c x + X j + k + l =4 c jkl x j p k q l + O (5) . The same steps as above lead to leading-order to the two-parameter normal form f ( x, p, q ) = c xp + ( c + c p + c q ) x . (7)Again, this shows that a second parameter can genericallychange second-order to first-order phase transitions; cf. Fig. 2.We now give three examples of complex systems where ageneralization leads to a change from a first-order to a second-order phase transition. We explicitly relate each of the exam-ples to the abstract framework above. Transitions in adaptive epidemics.
We consider the adap-tive epidemic model by Gross et al. [6], which is microscop-ically modeled as a Markov chain on networks with nodesbeing in two states, either susceptible S or infected I . Infec-tions take place at rate ρ , recovery at rate r (which we set to r = 1 without loss of generality here), and adaptive re-wiringof an SI -link to an SS -link at rate q . Direct numerical sim-ulations show that the bifurcation at the epidemic threshold ρ = ρ c is a second-order transition if q = 0 . It becomes afirst-order transition if q is increased sufficiently, i.e., the net-work becomes more strongly adaptive. Based upon our con-siderations above, it is natural to expect that allowing for gen-eral network topologies via re-wiring is a sufficiently genericbreaking mechanism to allow the second-to-first order changevia the parameter q . In fact, this is what is verified implic-itly in [6] by using a moment-closure expansions [10] of thenetwork dynamics. The following moment-closed ODEs de-scribe the dynamics for large networks I ′ = ρ ( µ − l II − l SS ) − I,l ′II = ρ ( µ − l II − l SS ) (cid:18) µ − l II − l SS − I + 1 (cid:19) − l II ,l ′SS = (1 + q )( µ − l II − l SS ) − ρ ( µ − l II − l SS ) l SS − I , where I and l II , l SS are a normalized infected density andtwo similarly normalized link densities respectively [6]; notethat conservation laws allow for the elimination of S and l SI .We fix µ arising from a connectivity assumption [6] of thenetwork to µ = 20 . This is a standard assumption [11], aswe only want to demonstrate the principal effect of adding re-wiring via q . It can be checked, see [6, 11], that a first-ordertransition is possible upon varying q .We now formally show that the change of criticality is aspecial case of our more general results above. One checksthat there always exists the invariant trivial branch of steadystates { I = 0 , l II = 0 , l SS = µ } . The epidemic thresholdbifurcation point is given by ρ c = 1 + qµ = 1 + q . Now we employ a lengthy, yet very direct and general, centermanifold calculation to find the normal form, which we out-line here. First, we shift coordinates I = X , l II = X , l SS = X + 10 , ρ = p + ρ c , to obtain a vector field X ′ = F ( X, p, q ) . Then we transform the linear part A = D X F (0 , , q ) into Jordan canonical form M − AM = − ( − q − for a transformation matrix X = M ˜ x that can be calculatedfrom the eigenvectors of A . We augment the new ODEs ˜ x ′ = M − F ( M ˜ x, p, q ) by p ′ = 0 and q ′ = 0 to cal-culate a the three-dimensional center manifold { (˜ x , ˜ x ) = h (˜ x , p, q ) } as there are three zero eigenvalues. The manifoldis parametrized over the center directions (˜ x , p, q ) . Usingthe invariance equation [8] and a quadratic ansatz for h , oneobtains after equating coefficients ˜ x = h (˜ x , p, q ) = − x , ˜ x = h (˜ x , p, q ) = 3364206763 ˜ x + 2801681 ˜ x p. Plugging this back into the equation for ˜ x ′ and writing x :=˜ x gives the flow on the center manifold to leading order as x ′ = 800 q + 41 xp + 80 (cid:16) q +1) q +41 − ( q + 1) (cid:17) q + 41 x + · · · Now one easily checks from the coefficients of xp and x ,as above, that the parameter q indeed yields a change in thecriticality from a second-order to a first-order transition at q =21 / . Hence, from this perspective we can clearly see that achange in criticality is not surprising: The re-wiring q appearsin the reduced center manifold as a sufficiently generic secondunfolding parameter as in the universal route described above. Synchronization in phase oscillator networks.
The Ku-ramoto model [12] describes the evolution of a network of N phase oscillators, where the state of oscillator k ∈ { , . . . , N } is given by the phase θ k ∈ R / (2 π Z ) . Kuramoto oscilla-tors interact additively, so a natural generalization is to con-sider the effect of non-additive interactions [13–15] since theyarise naturally in phase reductions of coupled nonlinear oscil-lators [16]. For example, Skardal and Arenas [4] consideredthe synchronization transition in such a variation of the Ku-ramoto model with triplet interactions. Specifically, the phaseof oscillator k evolves according to θ ′ k = ω k + K N N X j =1 sin( θ j − θ k )+ K N N X j,l =1 sin(2 θ l − θ j − θ k ) , (8)with intrinsic frequencies ω k sampled from a Lorentzian dis-tribution with mean and width [17]. The parameter K determines the strength of the additive interactions and K the strength of the triplet interactions; for K = 0 we recoverthe classical Kuramoto model.A sufficiently large triplet coupling strength K can nowchange the nature of the synchronization transition. Writei := √− and let Z = R e i φ = N P Nj =1 e i θ j denote the Ku-ramoto order parameter. In the mean-field limit of N → ∞ oscillators, the Ott–Antonsen reduction [18] of (8) yields theeffective dynamics R ′ = (cid:18) K − (cid:19) R + (cid:18) K − K (cid:19) R − K R , which describe the evolution of the system; see also [19]. Set p = K − , q = K and x = R to read off the normalform expansion (7). There is a critical transition at the bi-furcation p = 0 , that is K = 2 for any K . For the Ku-ramoto model K = 0 , the bifurcation is always supercritical(second-order). However, for K > the synchronizationtransition becomes subcritical and discontinuous. Hence, adiscontinuous synchronization transition in phase oscillatorswith higher-order interactions is another special case of theuniversal route described above. Discontinuous percolation transitions.
Change from acontinuous to a discontinuous transition have also been ob-served in percolation problems. Consider the q -state Pottsmodel on a Bethe-lattice with coordination number 3 [20]. For q = 2 this gives the Ising model. Now suppose that the bondsare occupied occupied independently with homogeneous den-sity ˆ p ∈ [0 , . Evaluating the percolation probabilities re-cursively [21], one obtains the percolation probability for thelattice as a fixed point of the iteration P n +1 = 2ˆ pP n + ( q − p P n p ( q − P n =: H ( P n ) . (9)The percolation transition of the fixed point P ∗ = 0 of H happens at the critical bond density ˆ p = ; whether this tran-sition is continuous or discontinuous depends on the numberof states q of the Potts model [21].The change of criticality of the percolation transition can beunderstood within the general framework introduced above.Set p := ˆ p − and consider the ODE x ′ = f ( x, p, q ) := H ( x ) − x (10)obtained by seeing (9) as a difference equation. By defini-tion, the fixed points of H in (9) correspond to equilibria of f in (10). Moreover, since ∂ x f = ∂ x H − and ∂ x H > ina neighborhood of ( x, p ) = 0 , linear stability of stationarystates coincides as well. Thus, the behavior of the percolationtransition of (9) is completely determined by the bifurcationsof the equilibrium x = 0 of (10) at p = 0 . A Taylor expansionof f ( x, p, q ) yields f ( x, p, q ) = 2 xp + 14 ( g ( p ) q − g ( p )) x + · · · with g ( p ) = 4 p +4 p +1 . Thus, the change of criticality of thepercolation transition at q = 2 corresponds to a change from asuper- to a subcritical transcritical bifurcation in the universalroute described above.How the the percolation probability changes with the bonddensity is also directly related to discontinuous transitions inthe expected maximal cluster size of a random graph. Specif-ically, random graphs with an underlying hierarchical self-similar structure allow to calculate the percolation probabil-ity through recursive relations [22] as in the Potts model dis-cussed above. By calculating the corresponding generatingfunctions [23], one can observe a discontinuous transition inthe expected size of the largest cluster. Discussion.
Our argument shows that from the perspectiveof bifurcation theory, one can expect a change from a con-tinuous transition to a discontinuous transition as additionaleffects are added to a classical model. Here, we gave threeexplicit network dynamics examples to illustrate this univer-sal route. Our formalism not only shows that a transition toan explosive change of system properties is not surprising butalso has the same underlying dynamical mechanism. In par-ticular, our framework links explosive percolation, explosivesynchronization, and explosive epidemic spreading explicitly.Many other variations are possible that fit into our frame-work. First, the effect of numerous generalizations of the Ku-ramoto model on the synchronization transition have been ex-plored in the literature. This includes varying the properties ofthe intrinsic frequencies [2, 3] or generalized coupling struc-tures that encode higher-order effects [4, 5]. Second, we an-ticipate our theory to be relevant in neural networks. For ex-ample, networks of quadratic-integrate-and-fire neurons canbe described by low-dimensional equations using a reductionclosely related to the Ott–Antonsen approach [19, 24]. Theseequations show transcritical bifurcation in a limiting case [25]that could shed light on the emergence of discontinuous tran-sitions between low- and high firing dynamics [26]. Finally,we expect our theory to apply also in further physical systems,for example, chemical reaction networks, where the same typeof mechanism is bound to be relevant.
Acknowledgments.
CB and CK gratefully acknowledgethe support of the Institute for Advanced Study at the Tech-nical University of Munich through a Hans Fischer Fellow-ship awarded to CB that made this work possible. CK ac-knowledges support via a Lichtenberg Professorship as wellas support via the TiPES project funded the European UnionsHorizon 2020 research and innovation programme under grantagreement No. 820970. [1] S. Boccaletti, J. Almendral, S. Guan, I. Leyva,Z. Liu, I. Sendi˜na-Nadal, Z. Wang, and Y. Zou,Phys Rep , 1 (2016); R. M. D’Souza, J. G´omez-Garde˜nes,J. Nagler, and A. Arenas, Adv in Phys , 123 (2019).[2] D. Paz´o and E. Montbri´o, Phys Rev E , 046215 (2009).[3] O. E. Omel’chenko and M. Wolfrum,Phys Rev Lett , 164101 (2012).[4] P. S. Skardal and A. Arenas, arXiv:1909.08057 (2019).[5] A. P. Mill´an, J. J. Torres, and G. Bianconi,arXiv:1912.04405 (2019).[6] T. Gross, C. D. D’Lima, and B. Blasius, Phys Rev Lett ,(208701) (2006).[7] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nat Commun ,2485 (2019).[8] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dy-namical Systems, and Bifurcations of Vector Fields (Springer,New York, NY, 1983).[9] H. Kielhoefer,
Bifurcation Theory: An Introduction with Appli-cations to PDEs (Springer, 2004).[10] C. Kuehn, in
Control of Self-Organizing Nonlinear Systems ,edited by E. Sch¨oll, S. Klapp, and P. H¨ovel (Springer, 2016) pp. 253–271.[11] C. Kuehn, J Nonlinear Sci , 457 (2013).[12] S. H. Strogatz, Physica D , 1 (2000); J. Acebr´on,L. Bonilla, C. P´erez Vicente, F. Ritort, and R. Spigler,Rev Mod Phys , 137 (2005).[13] M. Rosenblum and A. Pikovsky,Phys Rev Lett , 064101 (2007).[14] T. Tanaka and T. Aoyagi, Phys Rev Lett , 224101 (2011).[15] C. Bick, Phys Rev E , 050201(R) (2018).[16] P. Ashwin and A. Rodrigues, Physica D , 14 (2016); I. Le´onand D. Paz´o, Phys Rev E , 012211 (2019).[17] This assumption can be made without loss of generality by scal-ing time appropriately and going in a suitable co-rotating refer-ence frame.[18] E. Ott and T. M. Antonsen, Chaos , 037113 (2008). [19] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens,arXiv:1902.05307 (2020).[20] F. Y. Wu, Rev Mod Phys , 235 (1982); F. Peruggi, F. di Lib-erto, and G. Monroy, J Phys A-Math Gen , 811 (1983).[21] J. T. Chayes, L. Chayes, J. P. Sethna, and D. J. Thouless,Commun Math Phys , 41 (1986).[22] S. Boettcher, J. L. Cook, and R. M. Ziff,Phys Rev E , 041115 (2009).[23] S. Boettcher, V. Singh, and R. M. Ziff,Nat Comm , 787 (2012).[24] D. Paz´o and E. Montbri´o, Phys Rev X , 011009 (2014).[25] D. Paz´o and E. Montbri´o, Phys Rev Lett , 238101 (2016).[26] H. Schmidt, D. Avitabile, E. Montbri´o, and A. Roxin,PLOS Comp Bio14