A universal route to explosive phenomena

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 (11, 12) 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 ∈ ℝ^m 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 ∈ ℝ^m.

Now, suppose that the transition point corresponds to a bifurcation point (13) upon parameter variation. Generically, we can assume that a single eigenvalue of the Jacobian D_xF(0, y) ∈ ℝ^{n × n} crosses the imaginary axis. First, by a translation in parameter space, we may assume that the main bifurcation parameter is p = y₁. Second, using center manifolds (13) or Lyapunov-Schmidt reduction (14), 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, (15, 16).

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 ℤ₂ 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., (17). 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) = ρX² + aX⁴, 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.

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 = y₂, 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_jkl 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 c_0kl = 0 for all k, l ∈ ℕ₀ = ℕ ∪ {0}. Second, because a single eigenvalue crosses at p = 0, we must have ∂_xf(0,0, q) = 0, where ∂_x denotes the partial derivative with respect to x. Hence, we have c_10l = 0 for all l ∈ ℕ₀. Third, because we assume a simple eigenvalue crosses transversally, we get ∂_xp(0,0, q) ≠ 0 entailing c₁₁₀ ≠ 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 c₁₀₂ = c_0kl = 0 from above and c_jkl ≠ 0 if j ≥ 1. The leading-order nonvanishing conditions are ∂_xxpf(0) ≠ 0 and ∂_xxqf(0) ≠ 0, and we note that c₁₁₁xpq is of higher order in comparison to c₁₁₀xp for the linear part in x because c₁₁₀ ≠ 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 ∂_xxf(0, p, q) upon varying q in an interval [−q₀, q₀] for some q₀ > 0 as long as c₂₀₁ ≠ 0, which we expect generically as it is the leading-order term involving the parameter q. Even if c₂₀₁ = 0, we can expand to higher order in q and may thereby eventually change the sign of ∂_xxf(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 ∂_xxf(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.

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.

Epidemic dynamics on complex networks has been a very active topic for several decades (18). 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 (19), 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. (9) 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 (9) using a moment-closure expansions (11) of the network dynamics. The dynamics for large networks are described by the moment-closed ODEswhere I, and l_ℐℐ, l_SS are a normalized infected density and two similarly normalized link densities, respectively (9); note that conservation laws allow for the elimination of S and l_Sℐ. We fix μ arising from a connectivity assumption (9) of the network to μ = 20. This is a standard assumption (20), as we only want to demonstrate the principal effect of adding rewiring via q. It can be checked [see (9, 20)] 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 $\{I=0,l_{\mathcal{I}\mathcal{I}}=0,l_{\mathcal{S}\mathcal{S}}=\frac{\mathrm{\mu }}{2}\}$. The epidemic threshold bifurcation point is given byWe now compute the normal form using a direct and general center manifold calculation [see (13) for a general reference], which we outline here: First, we shift coordinates, I = X₁, l_ℐℐ = X₂, l_SS = X₃ + 10, ρ = p + ρ_c, to obtain a vector field $\stackrel{·}{X}=F(X,p,q)$. Then, we transform the linear part A = D_XF(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 $\stackrel{·}{z}=M^{-1}F(\mathit{Mz},p,q)$ by $\stackrel{·}{p}=0$ and $\stackrel{·}{q}=0$ to calculate the three-dimensional center manifold {(z₂, z₃) = h(z₁, p, q)} as there are three zero eigenvalues. The manifold is parameterized over the center directions (z₁, p, q). Using the invariance equation (13) and a quadratic ansatz for h, one obtains after equating coefficientsSubstituting this back into the equation for ${\stackrel{·}{z}}_1$ and writing x ≔ z₁ give the flow on the center manifold to leading order asThe coefficients of xp and x² 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 (21) or chemical reaction processes (22); 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.

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 (23, 24). For a network of N Kuramoto oscillators, the state of oscillator k ∈ {1, …, N} is given by the phase variable θ_k ∈ ℝ/(2πℤ). The phases evolve according tofor k = 1, …, N, where K₂ is the coupling strength between oscillators and the intrinsic frequencies ω_k sampled from a unimodal distribution. Write $i≔\sqrt{-1}$. The (complex-valued) Kuramoto order parameter $Z=R{\mathrm{e}}^{i\mathrm{\varphi }}=\frac{1}{N}{\sum }_{j=1}^N{e}^{i{\mathrm{\theta }}_j}$ 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 K₂ between oscillators is increased (25).

Kuramoto oscillators interact pairwise, so a natural generalization is to consider the additional effect of nonpairwise interactions (26–28) because they arise naturally in phase reductions of coupled nonlinear oscillators (29, 30). Skardal and Arenas (6) 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 ω_k 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 K₂ determines the strength of the additive interactions and K₃ determines the strength of the triplet interactions. While K₃ = 0 yields the classical Kuramoto model (Eq. 10) with a continuous synchronization transition, for sufficiently large K₃, 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 $\overline{w}$ 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 (12, 31): In the limit, the network dynamics are (exactly) described by the ODEof the order parameter Z; the derivation is given explicitly in (6). Substituting polar coordinates Z = Re^iϕ, 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 $p=\frac{K_2}{2}-1$, $q=\frac{K_3}{2}$, and x = R directly yields the normal form expansion (Eq. 9) of the pitchfork bifurcation. The bifurcation of R = 0 happens at p = 0 (K₂ = 2) for any K₃. For the Kuramoto model K₃ = 0, the pitchfork bifurcation is always supercritical (second order). However, for K₃ > 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 (32) are likely to provide further examples of the universal route. For example, Zhang et al. (8) 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 (33). 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 (4, 5) or generalized coupling structures that encode higher-order effects (6, 7).

The flow of fluids through porous media is an example of a percolation problem (34). In bond percolation, each link of an infinite lattice of nodes is occupied with probability $\widehat{p}$, and connected nodes form clusters of connected nodes. The percolation probability $P*(\widehat{p})$ describes whether a given node is part of an infinite cluster. Varying $\widehat{p}$, the system undergoes a critical transition in terms of $P*(\widehat{p})$—the percolation transition—as large-scale clusters emerge. By a classical result by Fortuin and Kasteleyn (35), these percolation problems can also be understood in terms of the Potts model (36) 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 (37). 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 $\widehat{p}\in [0,1]$. Evaluating the percolation probabilities recursively as a hierarchy of finite latices with percolation probability P_n, one obtains the percolation probability for the (infinite) lattice as a fixed point of the iterationas calculated in (38). The percolation transition of the fixed point P* = 0 of H happens at the critical bond density $\widehat{p}=\frac{1}{2}$; whether this transition in varying $\widehat{p}$ is continuous or discontinuous depends on the number of states q of the Potts model (38): 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 $p≔\widehat{p}-\frac{1}{2}$ 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 ∂_xf = ∂_xH − 1 and ∂_xH > 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) = 4p² + 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 (39) as in the Potts model discussed above. By calculating the corresponding generating functions (40), 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 (33) by varying an additional parameter (41).