A universal route to explosive phenomena.

Critical transitions are observed in many complex systems. This includes the onset of synchronization in a network of coupled oscillators or the emergence of an epidemic state within a population. "Explosive" first-order transitions have caught particular attention in a variety of systems when classical models are generalized by incorporating additional effects. Here, we give a mathematical argument that the emergence of these first-order transitions is not surprising but rather a universally expected effect: Varying a classical model along a generic two-parameter family must lead to a change of the criticality. To illustrate our framework, we give three explicit examples of the effect in distinct physical systems: a model of adaptive epidemic dynamics, for a generalization of the Kuramoto model, and for a percolation transition.

INTRODUCTION

Many complex nonlinear systems—ranging from epidemic spreading and synchronization of coupled oscillators to percolation on a network—undergo critical order-disorder transitions as a system parameter is varied. As in classical statistical mechanics (), these transitions can be continuous (second order) or discontinuous (first order) at the transition point. Discontinuous first-order transitions have attracted particular attention (, ) as they can lead to an “explosive” change of system properties. For a wide variety of complex systems, it has been observed that a variation of the model via additional features leads to the change from a continuous second-order to a discontinuous first-order critical transition. For example, the classical Kuramoto model shows a continuous synchronization transition. However, varying the distribution of intrinsic frequencies (, ), generalizing the network to simplicial or higher-order coupling (, ), or adding a dynamical rule that suppresses synchronization () all allow for discontinuous synchronization transitions. Similarly, adding adaptation () or higher-order coupling structures () to models of epidemic spreading can induce a discontinuous transition to the epidemic state. These and other examples follow the same paradigm: First, an additional effect is added to a classical model. Second, variation of a parameter associated with the new effect turns a previously second-order transition into a first-order transition.

Here, we give a mathematical argument that a change from a continuous to a discontinuous critical transition is not surprising in nonlinear dynamical systems but a generically/universally expected effect if additional parameters are varied. Specifically, we show that any typical model variation along a two-parameter family of any classical model with a second-order transition must lead to a change of the criticality to first order. First, this result shows that a change of the criticality of transition points in different complex systems has a common dynamical origin; we illustrate this in explicit examples involving adaptive epidemic dynamics, synchronization in the Kuramoto model with nonadditive higher-order interactions, and a model from percolation theory. Second, this insight can be useful to identify system perturbations to induce or prevent the emergence of abrupt critical transitions in a variety of physical systems. Third, it highlights that means to delay the onset of a critical transition—such as adding a suppressive rule—can change the nature of a critical transition from being continuous to being discontinuous.

RESULTS

A universal mechanism that modulates transitions

Normal forms for critical transitions

Here, we take a dynamical approach at critical transitions where the macroscopic dynamics change qualitatively; we relate the Landau’s classical approach to phase transitions below. Consider a high-dimensional dynamical system close to a critical transition point. Mean-field approximations are a commonly used tool to simplify such a system to a low-dimensional description in terms of mean-field variables. These approximations can be obtained through moment closure and other approaches; see (, ) and the examples below. Consequently, we assume that the mean-field or continuum limit dynamics of the physical system in question near the transition point are given by an ordinary differential equation (ODE)where x = x(t) ∈ ℝ denotes the mean-field variables and y ∈ ℝ denotes the model parameters. While the state undergoing a transition in the full system can be quite general—that is, it does not necessarily have to be a stationary solution of the microscopic dynamics—we assume that it corresponds to an equilibrium point for the macroscopic Eq. 1: Hence, suppose that x* = x*(y) is the corresponding smooth family of equilibrium points parameterized by y that exist for all parameters. Using a translation, we may shift the equilibria and assume that x* = 0 ≔ (0,0, …,0)⊤ ∈ ℝ is the equilibrium point, i.e., F(0, y) = 0 for all y ∈ ℝ.

Now, suppose that the transition point corresponds to a bifurcation point () upon parameter variation. Generically, we can assume that a single eigenvalue of the Jacobian D(0, y) ∈ ℝ crosses the imaginary axis. First, by a translation in parameter space, we may assume that the main bifurcation parameter is p = y1. Second, using center manifolds () or Lyapunov-Schmidt reduction (), the dynamics of the full system (Eq. 1) are locally given by the one-dimensional dynamics on the center manifoldwith a family of equilibria x*(p) = 0, that isthat bifurcate at p = 0. Note that center manifold reductions do not necessarily require a finite-dimensional approximation (Eq. 1) because it can also be directly applied to infinite-dimensional systems, such as partial differential equations or systems with delay; see, for example, (, ).

Now, it is well known in bifurcation theory that the two typical bifurcation points encountered in applications are the transcritical bifurcation with local normal formwhere a = ± 1 determines whether the bifurcation/transition upon varying p is second order (x ≥ 0, a = − 1) or first order (x ≥ 0, a = +1). Similarly, if there is an equivariance given by a ℤ2 reflection symmetry in the model via f(x, p) = − f( − x, p), then there cannot be any terms of even power in x in the Taylor expansion, and the generic transition is a pitchfork bifurcationThe pitchfork is second order if it is supercritical and a = −1, while it is first order if it is subcritical and a = +1.This dynamical perspective directly relates to Landau’s classical theory of phase transitions; cf., (). Briefly, this approach relies on defining an order parameter X and constructing a functional G(X) whose minima X* are the values the free energy of the system takes. For the Ising model of interaction spins, the order parameter X is the average of all spins, and the free energy functional is G(X) = ρX2 + aX4, where ρ is a shifted temperature and a ≠ 0 a parameter. By differentiation, we have that the minima X* that determine the free energy satisfywhich corresponds to equilibria of the normal form of the pitchfork bifurcation (Eq. 5). Hence, Landau’s phase transitions can be analyzed by looking at the branching behavior of equilibria (Eq. 3) in a corresponding dynamical model.

Change of criticality for a generic variation of an additional parameter

We now consider the variation of an additional parameter that takes into account the additional effect for each model as indicated above. Take a generic single eigenvalue crossing and with the phase space {x ≥ 0}. As the model is varied, the persistence of a single eigenvalue crossing is generic within one-parameter families of the vector field f. Hence, we may assume (without loss of generality) that the single eigenvalue crosses at p = 0. Furthermore, if we vary the model, then at least one additional free parameter, say q = y2, generically appears.

To understand how varying the additional parameter affects the equilibria, we now expand the vector field (Eq. 2) in x and p as well as q. A Taylor expansion at the bifurcation point yieldswhere O(M + 1) denotes terms of order M + 1. The coefficients c are constrained by the conditions imposed by the equilibrium and the bifurcation scenario we consider. First, the existence of a trivial branch of equilibria, f(0, p, q) = 0, implies c0 = 0 for all k, l ∈ ℕ0 = ℕ ∪ {0}. Second, because a single eigenvalue crosses at p = 0, we must have ∂(0,0, q) = 0, where ∂ denotes the partial derivative with respect to x. Hence, we have c10 = 0 for all l ∈ ℕ0. Third, because we assume a simple eigenvalue crosses transversally, we get ∂(0,0, q) ≠ 0 entailing c110 ≠ 0. In summary, we have

A generic model variation with at least one additional free parameter now leads to a vector field f that allows for a change in criticality. With the bifurcation conditions incorporated into (Eq. 6), one may use bifurcation theory to unfold the singular point into a generic family. In particular, the next derivatives of the vector field at the bifurcation point should not vanish. Hence, for combinations with j + k + l = 3, we must have c102 = c0 = 0 from above and c ≠ 0 if j ≥ 1. The leading-order nonvanishing conditions are ∂(0) ≠ 0 and ∂(0) ≠ 0, and we note that c111xpq is of higher order in comparison to c110xp for the linear part in x because c110 ≠ 0. Truncating higher-order terms, this yields the lowest-order two-parameter unfolding normal formWe now apply a scaling (or geometric desingularization, or renormalization) with a small parameter ε > 0 through the transformation (x, p, q) ↦ (xεα, pεβ, qεγ). For the transcritical normal form (Eq. 7), we choose α = 1, β = −1, γ = −2 to obtain (upon a suitable time rescaling)Hence, one easily checks that there is a sign change of ∂(0, p, q) upon varying q in an interval [−q0, q0] for some q0 > 0 as long as c201 ≠ 0, which we expect generically as it is the leading-order term involving the parameter q. Even if c201 = 0, we can expand to higher order in q and may thereby eventually change the sign of ∂(0, p, q). So only certain situations, e.g., the presence of symmetries or nongeneric smooth functions, could lead to the preservation of the sign for all q ∈ ℝ: A function without any dependence on the second parameter q may be the most extreme case of nongenericity, but symmetries can also force specific Taylor coefficients to vanish. Once the sign of ∂(0, p, q) changes, this implies that, generically, the second parameter is able to change the transition from second order to first order or vice versa. Of course, from the viewpoint of the geometry of the bifurcation diagram, this is quite intuitive, as shown in Fig. 1, that a second generic parameter may change criticality.

Fig. 1

Sketch for the variation of a transcritical bifurcation for the phase space {x ≥ 0} and parameters (p, q) with primary parameter p and second generic unfolding parameter q.

Dashed lines indicate instability of the equilibrium, and solid lines indicate stability. The gray cases are first-order (subcritical) transitions, while the black diagrams are second-order (supercritical) transitions.

The situation for the pitchfork works very similarly except that an additional symmetry f(x, p, q) = − f( − x, p, q) has to be respected. This further constrains the coefficients of the Taylor expansion. Note that if this symmetry is broken, then we are in the transcritical case if there is still a trivial branch for all values of the parameters. Hence, we now assume that the symmetry holds. Taylor expansion as above gives for a bifurcation point with a single eigenvalue crossingThe same steps as above lead to leading order to the two-parameter normal formAgain, this shows that a second parameter can generically change a second-order to a first-order transition; cf., Fig. 2.

Fig. 2

Sketch for the variation of a pitchfork bifurcation for the phase space {x ≥ 0} and parameters (p, q) with primary parameter p and second generic unfolding parameter q.

Dashed lines indicate instability of the equilibrium, and solid lines indicate stability. The gray cases are first-order (subcritical) transitions, while the black diagrams are second-order (supercritical) transitions.

Discontinuous critical transitions from a universal perspective

We now give three explicit examples of complex nonlinear systems where a generalization leads to a change from a first-order to a second-order critical transition. While seemingly distinct and from a variety of contexts, we show that the transitions are related through the abstract framework above.

Transitions in adaptive epidemics

Epidemic dynamics on complex networks has been a very active topic for several decades (). Classical susceptible-infected-susceptible (SIS) models are microscopically modeled as a Markov chain on networks with nodes being in two states, either susceptible S or infected ℐ. Infections take place at rate ρ along network links and recovery at rate r (which we set to r = 1 without loss of generality here). The effect of the underlying network is crucial to understand such contact processes (), for example, to determine the fraction of infected individuals in the long run.

However, changing social contacts affect the network, and thus, network adaptivity—dynamics of the network that interact with the dynamics on the network—affects epidemic dynamics. The paradigmatic and widely used adaptive epidemic model by Gross et al. () considers SIS dynamics with additional adaptive rewiring of an Sℐ link to an SS link at rate q. While this rule may suppress infection (as Sℐ links are removed), it can also create highly connected clusters of susceptible nodes (as SS links are added). Direct numerical simulations show that the bifurcation at the epidemic threshold ρ = ρc is a second-order transition if q = 0. It becomes a first-order transition if q is increased sufficiently, i.e., the network becomes more strongly adaptive. On the basis of our considerations above, it is natural to expect that allowing for general network topologies via rewiring is a sufficiently generic breaking mechanism to allow the second- to first-order change via the parameter q. This is what is verified implicitly in () using a moment-closure expansions () of the network dynamics. The dynamics for large networks are described by the moment-closed ODEswhere I, and l, lSS are a normalized infected density and two similarly normalized link densities, respectively (); note that conservation laws allow for the elimination of S and lSℐ. We fix μ arising from a connectivity assumption () of the network to μ = 20. This is a standard assumption (), as we only want to demonstrate the principal effect of adding rewiring via q. It can be checked [see (, )] that a first-order transition is possible upon varying q.

We now formally show that the change of criticality is a special case of our more general results above. One checks that there always exists the invariant trivial branch of steady states . The epidemic threshold bifurcation point is given byWe now compute the normal form using a direct and general center manifold calculation [see () for a general reference], which we outline here: First, we shift coordinates, I = X1, l = X2, lSS = X3 + 10, ρ = p + ρc, to obtain a vector field . Then, we transform the linear part A = D(0,0, q) into Jordan canonical formfor a transformation matrix X = Mz that can be calculated from the eigenvectors of A. We augment the new ODEs by and to calculate the three-dimensional center manifold {(z2, z3) = h(z1, p, q)} as there are three zero eigenvalues. The manifold is parameterized over the center directions (z1, p, q). Using the invariance equation () and a quadratic ansatz for h, one obtains after equating coefficientsSubstituting this back into the equation for and writing x ≔ z1 give the flow on the center manifold to leading order asThe coefficients of xp and x2 now show that the parameter q indeed yields a change in the criticality from a second-order to a first-order transition at q = 21/19. Hence, from this perspective, we can clearly see that a change in criticality occurs as expected: The rewiring q appears in the reduced center manifold as a sufficiently generic second unfolding parameter as in the universal route described above.

Of course, there are several other possible variations of this theme for concrete models of contact process character that use a secondary unfolding parameter to allow for a change criticality. Examples are epidemic models with additional node states () or chemical reaction processes (); the change of criticality corresponds to a “tricritical point.” However, our theoretical results clearly show that one should be able to go far beyond classical contact processes or reaction processes, as only needs the existence of a (local) dimension reduction or mean-field technique in combination with a suitable two-parameter bifurcation unfolding.

Synchronization in phase oscillator networks

The Kuramoto model has been instrumental to understand the emergence of synchrony in coupled oscillator networks, ranging from synchronization of flashing fireflies to emergent neural synchrony (, ). For a network of N Kuramoto oscillators, the state of oscillator k ∈ {1, …, N} is given by the phase variable θ ∈ ℝ/(2πℤ). The phases evolve according tofor k = 1, …, N, where K2 is the coupling strength between oscillators and the intrinsic frequencies ω sampled from a unimodal distribution. Write . The (complex-valued) Kuramoto order parameter describes the mean field of the oscillators. Specifically, its absolute value encodes the level of synchrony in the system: R = 0 if the phases of all oscillators are evenly distributed around the circle, and R = 1 if all oscillators are phase synchronized. This classical model exhibits—in analogy to phase transitions in statistical mechanics—a second-order transition from an incoherent state to a partially coherent state as the coupling strength K2 between oscillators is increased ().

Kuramoto oscillators interact pairwise, so a natural generalization is to consider the additional effect of nonpairwise interactions (–) because they arise naturally in phase reductions of coupled nonlinear oscillators (, ). Skardal and Arenas () showed that first-order transitions to synchrony arise in a variation of the Kuramoto model with nonadditive triplet interactions where the phase of oscillator k evolves according toand the intrinsic frequencies ω are sampled from a Lorentzian distribution with mean 0 and width 1; the choice of parameter for the Lorentzian can be made without loss of generality by scaling time appropriately and going in a suitable corotating reference frame. The parameter K2 determines the strength of the additive interactions and K3 determines the strength of the triplet interactions. While K3 = 0 yields the classical Kuramoto model (Eq. 10) with a continuous synchronization transition, for sufficiently large K3, this transition can become discontinuous.

The change to a discontinuous synchronization transition in phase oscillators with higher-order interactions can be understood in terms of the universal route described above. Let denote the complex conjugate of a complex number w. In the mean-field limit of N → ∞ oscillators, the network dynamics (Eq. 11) can be described using the Ott-Antonsen reduction (, ): In the limit, the network dynamics are (exactly) described by the ODEof the order parameter Z; the derivation is given explicitly in (). Substituting polar coordinates Z = Re, the dynamics of mean-field phase ϕ, and mean-field amplitude R decouple, yielding the effectively one-dimensional dynamicsThe nature of the synchronization transition is determined by the bifurcation of the equilibrium R = 0. Setting , , and x = R directly yields the normal form expansion (Eq. 9) of the pitchfork bifurcation. The bifurcation of R = 0 happens at p = 0 (K2 = 2) for any K3. For the Kuramoto model K3 = 0, the pitchfork bifurcation is always supercritical (second order). However, for K3 > 2, the synchronization transition becomes a subcritical first-order transition in line with our universal approach.

Many other variations of the Kuramoto model that change the nature of the synchronization transition () are likely to provide further examples of the universal route. For example, Zhang et al. () considered a modified Kuramoto model where the coupling depends on the intrinsic frequency. This yields coupling that suppresses oscillator synchronization as a dynamical analog to the Achlioptas rule for percolation in random networks (). The Kuramoto model with this additional feature exhibits a discontinuous synchronization transition and relates our theory with explosive synchronization. Other variations of the Kuramoto model that affect the synchronization transition include varying the properties of the intrinsic frequencies (, ) or generalized coupling structures that encode higher-order effects (, ).

Discontinuous percolation transitions

The flow of fluids through porous media is an example of a percolation problem (). In bond percolation, each link of an infinite lattice of nodes is occupied with probability , and connected nodes form clusters of connected nodes. The percolation probability describes whether a given node is part of an infinite cluster. Varying , the system undergoes a critical transition in terms of —the percolation transition—as large-scale clusters emerge. By a classical result by Fortuin and Kasteleyn (), these percolation problems can also be understood in terms of the Potts model () as a generalization of the Ising model. This correspondence allows one to relate the percolation transition and phase transitions in the Potts model.

Whether the percolation transition is continuous or discontinuous now depends on the system parameters. For the Potts model on a Bethe lattice, the percolation probability can be evaluated using recursive relations (). Consider the q-state Potts model on a Bethe lattice with coordination number 3. For q = 2, this gives the Ising model. Now, suppose that the bonds are occupied independently with homogeneous density . Evaluating the percolation probabilities recursively as a hierarchy of finite latices with percolation probability P, one obtains the percolation probability for the (infinite) lattice as a fixed point of the iterationas calculated in (). The percolation transition of the fixed point P* = 0 of H happens at the critical bond density ; whether this transition in varying is continuous or discontinuous depends on the number of states q of the Potts model (): If 1 ≤ q ≤ 2, then the transition is continuous, and if q > 2, then the transition is discontinuous.

This change of criticality of the percolation transition can be understood within the general dynamical framework introduced above. Set and consider the ODEobtained by seeing Eq. 12 as a difference equation. By definition, the fixed points of H in Eq. 12 correspond to equilibria of f in Eq. 13. Moreover, because ∂ = ∂ − 1 and ∂ > 0 in a neighborhood of (x, p) = 0, linear stability of stationary states coincides as well. Thus, the behavior of the percolation transition of Eq. 12 is completely determined by the bifurcations of the equilibrium x* = 0 of Eq. 13 at p = 0. A Taylor expansion of f(x, p, q) yieldswith g(p) = 4p2 + 4p + 1. Thus, the change of criticality of the percolation transition at q = 2 corresponds to a change from a supercritical to a subcritical transcritical bifurcation in the universal route described above.

How the percolation probability changes with the bond density is also directly related to discontinuous transitions in the expected maximal cluster size of a random graph. Specifically, random graphs with an underlying hierarchical self-similar structure allow one to calculate the percolation probability through recursive relations () as in the Potts model discussed above. By calculating the corresponding generating functions (), one can observe a discontinuous transition in the expected size of the largest cluster. Last, the critical transition in the q-state Potts model has also inspired a rule that links classical bond percolation and a type of Achlioptas process () by varying an additional parameter ().

DISCUSSION

Our argument shows that from the perspective of bifurcation theory, one can expect a change from a continuous transition to a discontinuous transition as additional effects are added to a classical model. Here, we gave three explicit physical examples to illustrate this universal route. Our formalism shows that a transition to an abrupt change of system properties is not only expected but also has the same underlying dynamical mechanism. In particular, our framework links the emergence of abrupt critical transitions for percolation, synchronization, and epidemic spreading explicitly through a common dynamical framework.

Many other variations—beyond the ones already discussed in the context of the explicit examples—are possible that fit into our framework. First, our mechanism also yields a natural explanation for first-order synchronization transitions observed in coupled nonlinear oscillators (), because higher-order effects affect the dynamics beyond the weak coupling limit (). Second, we anticipate our theory to be relevant in neural networks. For example, networks of quadratic-integrate-and-fire neurons can be described by low-dimensional equations using a reduction closely related to the Ott-Antonsen approach (, ). These equations show transcritical bifurcation in a limiting case () that could shed light on the emergence of discontinuous transitions between low and high firing dynamics (). Last, we expect our theory to apply also in further physical systems, for example, laser dynamics (), flocking (), or chemical reaction networks, where the same type of mechanism is bound to be relevant.