Explosive Phenomena in Complex Networks

TL;DR

This paper reviews explosive phenomena in complex networks, focusing on abrupt transitions like explosive percolation and synchronization, their models, properties, and implications across various real-world systems.

Contribution

It synthesizes and classifies existing models of explosive phenomena, highlighting their mechanisms, properties, and applications in complex networks.

Findings

Explosive percolation can be continuous with atypical finite-size scaling.

Models show that explosive synchronization relates to network structure and frequency distribution.

Explosive phenomena are linked to real-world events like epidemics and electric breakdowns.

Abstract

The emergence of large-scale connectivity and synchronization are crucial to the structure, function and failure of many complex socio-technical networks. Thus, there is great interest in analyzing phase transitions to large-scale connectivity and to global synchronization, including how to enhance or delay the onset. These phenomena are traditionally studied as second-order phase transitions where, at the critical threshold, the order parameter increases rapidly but continuously. In 2009, an extremely abrupt transition was found for a network growth process where links compete for addition in attempt to delay percolation. This observation of "explosive percolation" was ultimately revealed to be a continuous transition in the thermodynamic limit, yet with very atypical finite-size scaling, and it started a surge of work on explosive phenomena and their consequences. Many related models are now shown to yield discontinuous percolation transitions and even hybrid transitions. Explosive percolation enables many other features such as multiple giant components, modular structures, discrete scale invariance and non-self-averaging, relating to properties found in many real phenomena such as explosive epidemics, electric breakdowns and the emergence of molecular life. Models of explosive synchronization provide an analytic framework for the dynamics of abrupt transitions and reveal the interplay between the distribution in natural frequencies and the network structure, with applications ranging from epileptic seizures to waking from anesthesia. Here we review the vast literature on explosive phenomena and synthesize the fundamental connections between models and survey the application areas. We attempt to classify explosive phenomena based on underlying mechanisms and to provide a coherent overview and perspective for future research to address the many vital questions that remained unanswered.

Tables

Table 1. Table 1: Seminal models of discontinuous percolation transitions which pre-date the work on EP.

Model name	Process considered	First reference
Bootstrap percolation	Pruning of sites with too few neighbors	Chalupa et al. (1979)
$k$-core percolation	Pruning of insufficiently connected nodes	Seidman (1983); Bollobás (1984)
Rigidity percolation	Dynamics of rigid and floppy regions in glassy structures	Thorpe (1985)
Jamming percolation	Jamming transition in sphere packing	O’Hern et al. (2002)
Generalized contagion	SIR threshold dynamics	Dodds and Watts (2004)
Tricritical dynamic percolation	4-species spreading dynamics	Janssen et al. (2004)

Table 2. Table 2: The BFW algorithm . At each discrete time step u 𝑢 u , the randomly selected edge e u subscript 𝑒 𝑢 e_{u} is examined via this algorithm, where u 𝑢 u denotes the total number of edges sampled, A 𝐴 A the set of edges which have been accepted thus far (initially A = ∅ 𝐴 A=\emptyset ), and T = | A | 𝑇 𝐴 T=|A| the number of accepted edges.

Set $l=$ maximum size component in $A\cup \left\{e_u\right\}$
if $(l\le k)\{$
$A\leftarrow A\cup \left\{e_u\right\}$
$u\leftarrow u+1.$ (Get next edge.)}
else if $(T/u\ge g\left(k\right))\{u\leftarrow u+1$. (Get next edge.)}
else { $k\leftarrow k+1$. Then repeat this block.}

Table 3. Table 3: Key distinctions between EP and ES.

Explosive Percolation	Explosive Synchronization
Kinetic process	Dynamical process
Cluster nodes are linked	Cluster nodes are in sync
Typically non-local cluster merging	Local node-to-node-based cluster merging
Sufficiently homogeneous initial cluster size distribution	Heterogeneous node frequency distribution (sufficient condition)
Gap between small-sized and large-sized clusters (powder keg)	Frequency gap condition (Sec: III.6)

Equations

$$\chi =\frac{1}{N}\sum _j{j}^2{N}_j$$

$$\chi \sim |t-t_c{|}^{-\gamma }.$$

$$\frac{d{n}_k}{dt}=\sum _{i+j=k}{K}_{ij}{n}_i{n}_j-2n_k\sum _j{K}_{kj}{n}_j,$$

$${\Delta }_N(0.5,0.5)\propto N^{2/3}.$$

$${\Delta }_N(0.5,0.5)/N\propto N^{-1/3}\to 0\mathrm{as}N\to \infty$$

$$\Delta C_{\max }\equiv \mathrm{ma}{\mathrm{x}}_{\{\mathrm{T}\}}\left[|C(T+1)|-|C(T)|\right]/N.$$

$$W=\frac{1}{N}F\left(t(N^{\alpha }),N^{1-\beta }\right)$$

$$P(s,t)=\frac{sn(s)}{⟨s⟩},$$

$$P_{p=p_c,N}\left(\frac{s_{\max }}{N}\right)\sim N^{\eta }f\left(\frac{s_{\max }}{N}\cdot N^{\eta }\right).$$

$$\min \left\{1,\exp \left[-\alpha {\left(\frac{s-⟨s⟩}{⟨s⟩}\right)}^2\right]\right\},$$

$$g(k)=\alpha +(2k{)}^{-\beta }$$

$$K_{ij}=(ij{)}^{\omega }/(\sum _s{s}^{\omega }{n}_s{)}^2$$

$$P=\{\begin{array}{ll}0& \text{if }p<p_c\\ \Delta +f(p)& \text{if }p\ge p_c,\end{array}$$

$$P=\{\begin{array}{ll}0& \text{if }p<1\\ 1& \text{if }p=1.\end{array}$$

$$P=\{\begin{array}{ll}0& \text{if }p<p_c\\ \Delta +(p-p_c{)}^{\beta }& \text{if }p\ge p_c,\end{array}$$

$${\mathcal{R}}_v(|C|)=\frac{⟨|C|-⟨|C|⟩{⟩}^2}{⟨|C|{⟩}^2}$$

$$|C|/N\approx (t-t_c{)}^{\beta }=(1/N{)}^{\beta }.$$

$$N\le {\left(\frac{1}{\Delta }\right)}^{1/\beta }.$$

$${\dot{\theta }}_i={\omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j-{\theta }_i)(i=1,...,N),$$

$$r(t){\text{e}}^{\text{i}\varphi (t)}=\frac{1}{N}\sum _{j=1}^N{\text{e}}^{\text{i}{\theta }_j(t)},$$

$${\dot{\theta }}_i={\omega }_i+Kr\sin (\varphi -{\theta }_i)(i=1,...,N).$$

$$r=Kr{\int }_{\frac{-\pi }{2}}^{\frac{\pi }{2}}{\cos }^2\theta g(\omega )d\theta ,$$

$$K_c=\frac{2}{\pi g(0)}.$$

$${\dot{\theta }}_i={\omega }_i+\sum _j{\lambda }_{ij}{A}_{ij}\sin ({\theta }_j-{\theta }_i)(i=1,...,N),$$

$$r_i{\text{e}}^{\text{i}{\varphi }_i}=\sum _{j=1}^N{A}_{ij}⟨{\text{e}}^{\text{i}{\theta }_j}{⟩}_t,$$

$$r_i$$

$$r=\frac{r_i}{k_i}=\frac{1}{k_i}\left|\sum _{j=1}^N{A}_{ij}⟨{\text{e}}^{\text{i}{\theta }_j}{⟩}_t\right|.$$

$$\sum _j^N{k}_j=\lambda \sum _j^N{k}_j^2{\int }_{-1}^{1}g(x\lambda r{k}_j)\sqrt{1-x^2}dx.$$

$$\int kP(k)dk$$

$${\lambda }_c=\frac{2}{\pi g(0)}\frac{⟨k⟩}{⟨k^2⟩}=K_c\frac{⟨k⟩}{⟨k^2⟩}.$$

Full text

1. I Introduction

2. I.1 Percolation, synchronization, and phase transitions

3. I.2 Percolation and synchronization on complex networks

4. I.3 Overview of Review

5. II Explosive percolation (EP)

6. II.1 Percolation on random graphs

7. II.1.1 Static formulation

8. II.1.2 Kinetic formulation

9. II.1.3 Discontinuous percolation transitions

10. II.2 EP from competitive percolation

11. II.2.1 Achlioptas Process (AP)

12. II.2.2 AP with $m$-edge rules on random graphs

13. II.2.3 AP with $m$-edge rules on lattices

14. II.2.4 AP with $k$-vertex rules on random graphs

15. II.2.5 Thermodynamic formulations of EP

16. II.3 EP from alternate processes

17. II.3.1 Suppressing the largest component

18. II.3.2 Tuning edge rejection rate (The BFW model)

19. II.3.3 Cluster aggregation

20. II.3.4 Correlated processes

21. II.3.5 Hierarchical lattices

22. II.4 Categories of EP transitions

23. II.4.1 Type I: Anomalous continuous transition

24. II.4.2 Type II: Discontinuous phase transition

25. II.4.3 Type III: Discontinuity at the end of the process

26. II.4.4 Type IV: Single-step-continuous discontinuity

27. II.4.5 Type V: Hybrid phase transition

28. II.4.6 Type VI: Devil’s staircases and non-self-averaging

29. II.5 EP phenomena in real-world systems

30. II.5.1 Non-Self-Averaging and the emergence of molecular life

31. II.5.2 Finite size of real systems

32. II.5.3 Modular network structure in the real-world

33. II.5.4 Disordered materials

34. II.5.5 Explosive Immunization

35. II.5.6 Optimal Dismantling

36. II.5.7 Information spreading

37. II.6 Other interesting EP phenomena

38. II.6.1 Node arrival mitigates abruptness

39. II.6.2 APs with directed edges

40. II.6.3 Degree-based rules

41. II.6.4 Explosive Ising from an AP

42. II.6.5 Discrete Scale Invariance

43. II.7 General insights

44. II.7.1 Mechanisms and processes underlying EP

45. II.7.2 The role of order parameter

46. II.7.3 Reversible versus irreversible transitions

47. III Explosive Synchronization (ES)

48. III.1 Introduction to the Kuramoto model

49. III.2 The Kuramoto model in networks

50. III.3 Explosive Synchronization in networks

51. III.4 The correlated degree-frequency (CDF) framework

52. III.4.1 Analytical characterization

53. III.4.2 The star network

54. III.4.3 Effect of structural properties

55. III.4.4 More general CDF scenarios

56. III.4.5 Noisy degree-frequency correlations

57. III.4.6 CDF with time delay and frustration

58. III.5 The correlated frequency-coupling (CFC) framework

59. III.5.1 Independence of network properties

60. III.5.2 The Sakaguchi-Kuramoto model under the CFC framework

61. III.5.3 Loss of rotational invariance

62. III.5.4 Including frequency mismatch

63. III.5.5 Effective Synchronization centrality

64. III.6 The frequency gap conditioned (FGC) framework

65. III.7 The adaptive framework

66. III.8 Other paradigmatic dynamical models

67. III.8.1 Second-order Kuramoto model

68. III.8.2 Oscillators with amplitude dynamics

69. III.8.3 Neural models

70. III.8.4 Rossler chaotic oscillators

71. III.9 Experiments and real-world applications

72. III.9.1 Implementations in the lab

73. III.9.2 Spontaneous synchronization in real systems

74. III.10 ES in a nutshell

75. IV Connection between EP and ES

76. V Other models with explosive behavior

77. V.1 Explosive congestion in complex networks

78. V.2 Explosive epidemics in networks

79. V.3 Explosive phenomena on multilayer complex networks

80. V.3.1 Abrupt structural transitions

81. V.3.2 Abrupt percolation transitions

82. V.3.3 Abrupt synchronization transitions

83. VI Conclusions and Open Challenges

I Introduction

This review of explosive phenomena in complex networks is an attempt to provide a coherent overview of an area with vigorous activity and extensive literature dating back to 2009. We discuss the paradigms of Explosive Percolation (EP) and Explosive Synchronization (ES) in complex networks, surveying the broad range of phenomena observed in models that display these types of abrupt transitions. We identify that these models share a common underlying mechanism with a microscopic evolution dynamics that delays the formation of a macroscopic (connected or synchronized) component. We show that similar underlying cluster dynamics can be associated with many real systems displaying abrupt transitions. We also discuss examples of explosive phenomena in a variety of scenarios such the emergence of globalization from local individual decisions, epileptic seizures, explosive epidemics, and the emergence of molecular life.

I.1 Percolation, synchronization, and phase transitions

Percolation theory is used to analyze the extent of connectivity in an underlying random media such as a lattice or a random graph. It has been used extensively for many decades to model phenomena ranging from flow through porous media to connectivity and spreading on random graphs Stauffer and Aharony (1994); Sahimi (1994); Dorogovtsev et al. (2008). It was first formulated as a lattice model where each of the “bonds” between adjacent neighbors could be occupied with probability $p$ or unoccupied with probability $(1-p)$ Broadbent and Hammersley (1957). A cluster of adjacent occupied bonds is considered to be connected, allowing a fluid or an epidemic to spread throughout that cluster. Of great interest is the percolation phase transition that describes the sudden, but typically smooth, onset of large scale connectivity at a critical value $p_{c}$.

While percolation focuses on the structural properties of a system, synchronization constitutes a paradigm for analyzing its dynamical behavior. Synchronization is a fundamental phenomenon which is pervasive in systems ranging from coupled mechanical clocks to neuronal firing patterns in the brain. Although percolation and synchronization consider different aspects of a system, both phenomena are traditionally known to emerge via a second-order phase transition. In other words, starting from a disconnected system, long-range connectivity (i.e., percolation) emerges smoothly at a critical point when increasing internal system density. In turn, starting from a system of uncoupled identical oscillators, global synchronization emerges smoothly at a critical point when increasing the coupling between oscillators. In both cases, the global behavior (percolation or synchronization) arises as a consequence of increasing the strength of local interactions and are, in general, second-order phase transitions.

Thermal equilibrium phase transitions are typically classified as first-order or second-order based on how the order parameter changes at the critical point, with a discrete jump in the value for first-order, and a sharp but smooth change for second-order. Second-order transitions are additionally accompanied by a range of interrelated scaling behaviors, whereas first-order transitions typically display instead other phenomena like phase coexistence and hysteresis loops indicative of irreversibility. See, for instance, Stanley (1971), and the discussion therein.

In the scope of percolation, several examples of first-order transitions were known before the introduction of EP, such as $k$-core percolation Seidman (1983); Bollobás (1984). In the scope of synchronization, several works pointed out that the character of the phase transition may be first-order in specific setups that are unrelated to the underlying topology of connections, but instead are related to the particular distribution of natural frequencies, to random fields, or to inertial terms in the dynamical equations Kuramoto (1984); Bonilla et al. (1992); van Hemmen and Wreszinski (1993); Tanaka et al. 1997c ; Pazó (2005). Here our focus is on the more recently discovered phenomena of explosive transitions. These are abrupt transitions that arise as a consequence of the use of microscopic rules that, supported by the networked architecture, aim to hinder and postpone the formation of macroscopic components. In such a way, once the emergence of collective states becomes unavoidable, they emerge abruptly.

I.2 Percolation and synchronization on complex networks

With the increasing prevalence and reliance of modern society on a collection of networks, including social, transportation, biological, ecological, and communication networks, a coherent scientific study of the principles of networks started emerging in the late 1990s Watts and Strogatz (1998); Barabási and Albert (1999); Kleinberg et al. (1999). The paradigms of percolation and synchronization have provided theoretical underpinnings for analyzing the structure, function, resilience and robustness of network systems (see for instance Barrat et al. (2008); Newman (2010)). The classic starting points are the Erdős-Rényi model of a random network discussed in Sec. II.1 and the Kuramoto model of synchronization discussed in Sec. III.1.

In any network system, the extent of connectivity is a fundamental property that shapes what functions that network can support. For instance, ensuring large-scale connectivity is essential for transportation networks such as the world-wide airline network or for communication systems like the internet. Yet, when an infectious disease is spreading over a network, large-scale connectivity becomes a liability; a virus spreading on a well-connected social or computer network can reach enough nodes to cause an epidemic. Thus, in some contexts large-scale connectivity is highly desirable, yet in other contexts large-scale connectivity is quite hazardous.

The phenomena of “explosive percolation” was first discovered in an attempt to use small interventions to delay the onset of large-scale connectivity Achlioptas et al. (2009). The quest to understand the true nature of the resulting abrupt transition and the associated underlying mechanisms sparked a flurry of activity into alternative models, e.g., da Costa et al. (2010); Nagler et al. (2012); Riordan and Warnke (2012), thermodynamics of cluster growth, e.g., Cho et al. (2009), and explosive phenomena as a modeling paradigm, leading to the discovery of additional phenomena. A number of connections of EP to real-world systems have begun to emerge, e.g., Clusella et al. (2016).

In parallel, explosive synchronization in networks was first reported by Gómez-Gardeñes et al. 2011a when coupling the intrinsic frequencies of a dynamical system of phase oscillators with a topological characteristic of individual units, namely their number of local connections. An abrupt transition was observed in the order parameter that measures the global level of synchronization, for certain values of the coupling between units. This abrupt change is accompanied by an hysteresis cycle characteristic of first-order phase transitions, Sethna et al. (1993), and has been proven to exist as a bistability region in the phase space. These results have stimulated a large number of follow-on studies, e.g., Boccaletti et al. (2016), including establishing connections between ES and the phenomenology of several real systems, e.g. Kim et al. (2016); Wang et al. (2016). Similar to EP, in ES the emergence of large-scale collective behavior is delayed until it emerges dramatically.

I.3 Overview of Review

Our attempt here is to provide a coherent overview organizing the literature and identifying the most fruitful open directions for future work. The review is organized as follows. First in Sec. II we review Explosive Percolation phenomena. In Sec. III we review Explosive Synchronization phenomena. In Sec. IV we review the connections between EP and ES. In Sec. V we discuss other classes of dynamical processes on complex networks that show explosive transitions, including some recent works on multilayer and interdependent networks. Finally, in Sec. VI we provide conclusions and open challenges for the field.

Note that the scope of this review goes beyond being a summary and overview of the literature, and also provides a comprehensive analysis of the roots behind explosive behavior and the connections between the apparently different methods and models that display the associated phenomena. To this aim, in those parts devoted to EP and ES the reader will find brief introductions to the basics of percolation and synchronization phenomena on networks (due to the variety of phenomena covered in Sec. V, we do not devote such attention there). In addition, while the EP section describes in detail the statistical-physics based models and analysis of network growth leading to EP, the section devoted to ES is more focused on theoretical analysis of the equations underlying the nonlinear models leading to ES. The difference between both approaches is rooted in the differing natures of EP and ES and offers a complementary view of explosive transitions that we believe will pave the way for a deeper understanding of explosive phenomena.

II Explosive percolation (EP)

Explosive percolation (EP) describes the critical behavior of a general class of graph evolution models with microscopic dynamics that delay the growth of large components until large-scale connectivity inevitably emerges in a dramatic and abrupt manner, hence the term “explosive”. EP transitions were first introduced in Achlioptas et al. (2009), by a variant on standard percolation where a simple selection criteria is used to determine how edges (i.e., links) are added between a collection of $N$ originally isolated nodes (i.e., vertices). Much more is now understood about the broad range of models and behaviors that fit the EP paradigm. EP transitions have now been shown to exhibit an array of novel universality classes, and anomalous critical and supercritical behaviors. For instance, in the supercritical regime the order parameter is not necessarily a function of the control parameter; multiple discontinuous transitions including a “Devil’s staircase” of supercritical discontinuous transitions arbitrarily close to the initial percolation transition are possible; multiple giant components can coexist, which is not possible to classic percolation, and this feature also gives rise to modular structure.

Many underlying mechanisms that give rise to EP have now been identified including growth by overtaking, directly suppressing the largest component, preferentially keeping all clusters similar in size to the average, and using correlated percolation processes. It has been shown that EP processes can exhibit a discrete scale invariance predicting the location of the percolation transition. Many connections to real-world networks have been identified and EP is an emerging modeling paradigm. In this section we aim to synthesize these existing results and also the fundamental concepts and mechanisms underlying EP.

II.1 Percolation on random graphs

Here we briefly review the foundational models of percolation on random networks, including the paradigm of cluster aggregation and the mechanism of “multiplicative coalescence” that underlies the continuous emergence of a unique giant component. We refer to these models as classic random graph models of percolation.

II.1.1 Static formulation

Several related models defining the concept of a random graph Solomonoff and Rapoport (1951); Gilbert (1959); Erdős and Rényi (1959, 1960) were introduced in the 1950’s. They consider a collection of $N$ vertices with edges connecting vertices uniformly at random under two slightly different formulations. The first formulation considers that every possible edge has probability $p$ of being present (i.e., “occupied”), which is referred to as the $G(N,p)$ model Solomonoff and Rapoport (1951); Gilbert (1959). The second formulation considers that exactly $m$ edges are chosen to be occupied uniformly at random between the $N$ vertices, which is referred to as the $G(N,m)$ model Erdős and Rényi (1959, 1960). The two formulations are coincident in expectation, with the expected number of edges $\left<T\right>=m=pN(N-1)/2\approx pN^{2}/2$ (with the final approximation becoming exact as $N\rightarrow\infty$). Thus, in the large $N$ limit, the ratio of the expected number of edges to nodes is $t=m/N=pN/2$, which provides the control parameter to parameterize the process.

Nodes that are connected to one another following a path of occupied edges are considered to be in the same component (i.e., cluster) and we are interested in the distribution of component sizes for a given value of $t$. Of particular interest is the size of the largest connected component, $|C|$. This is an order parameter that displays a second order phase transition at $t_{c}=1/2$. For $t<t_{c}$, $|C|$ is of order logarithmic in $N$ and, for $t>t_{c}$, there is a unique largest component with size linear in $N$, with the value of $|C|$ changing rapidly but smoothly close to $t_{c}$ Bollobás (2001). In fact, for $t$ slightly greater than $t_{c}$ the size of the largest component is described by the function $|C|\sim(4t-2)N$  Bollobás (2001).

One can further characterize the nature of the random graph phase transition by analyzing associated response functions, such as what is referred to in the literature as the “susceptibility”, $\chi$, defined as the second moment of the component sizes, see e.g., Bollobás (2001); Grimmett (2010); Newman (2010). In more detail, this can be written

[TABLE]

where $N_{j}$ is the number of clusters that are of size $j$. Note $jN_{j}/N$ is the fraction of vertices in components of size $j$, so this second moment, $\chi$, is the expected size of the component to which an arbitrary vertex belongs. In the critical regime we see divergence of $\chi$ as expected for a standard second order phase transtion

[TABLE]

Classic random graph percolation obeys the mean-field exponents, with $\gamma=1$. For an explicit treatment, see for instance  Bollobás (2001); Saberi (2015). It should be noted that although this definition of susceptibility does show a divergence at the critical point, it does not obey the fluctuation-dissipation theorem that gives the connection between susceptibility and variances in statistical thermodynamics Bizhani et al. (2012).

II.1.2 Kinetic formulation

The model above can be recast as a kinetic process Ziff and Stell (1980); Ben-Naim and Krapivsky (2005) using the framework of cluster aggregation and the celebrated Smoluchowski equation von Smoluchowski (1916). Beginning from a collection of $N$ isolated nodes, edges are added, one at a time, uniformly at random to the graph in a discrete time process, with $T$ denoting the total number of edges added. The control parameter is once again the relative number of introduced edges $t=T/N$.

Since the process begins with a collection of $N$ isolated nodes, the first edge necessarily joins together two previously distinct components, forming a component of size two. Similarly, while the component sizes are sufficiently small, the likelihood that an edge chosen at random is internal to a component is negligible and with high probability each added edge joins together two previously disjoint components. With this assumption, one can then describe the random graph evolution process as a Smoluchowsky rate equation of cluster aggregation von Smoluchowski (1916). Here the likelihood that a random edge merges together a component of size $i$ with a component of size $j$ is described by a collision kernel $K_{ij}$. The evolution of the corresponding cluster density is described by the Smoluchowski rate equation as follows:

[TABLE]

where $n_{k}$ is the density of $k$-size clusters (the number of clusters of size $k$, divided by the total number of clusters). The particular collision kernel, $K_{ij}$, captures the adhesion properties and clusters’ diffusivity. The rate $K_{ij}n_{i}n_{j}$ reflects the probability that two clusters of sizes $i$ and $j$ merge per unit volume and unit time thus producing a cluster of size $k=i+j$. The negative term accounts for the case that a $k$-size cluster merges with any of the remaining clusters.

Cluster aggregation analysis is used in many studies discussed throughout this review. Modern perspective and rigorous details on the approach can be found in the recent textbook Krapivsky et al. (2010). The approach ignores intra-cluster link structure of the underlying network, but from the computational perspective, this allows for much faster simulations as only the distribution of cluster sizes needs to be propagated. In particular the celebrated Newman-Ziff algorithm based on union-find Newman and Ziff (2001) is pervasively used for efficient computation of percolation.

For the random graph model discussed thus far, two vertices are chosen uniformly at random and connected by an edge. The likelihood of choosing a vertex in a particular component of size $i$ is proportional to $i$ (the number of vertices in that component). Thus, the likelihood of merging two components is proportional to the product of their sizes and $K_{ij}\propto ij$, which is referred to as “multiplicative coalescence” Aldous (1999). This leads to a phase transition in the size of the largest component (called gelation) that is mathematically equivalent to static formulations of percolation on a random network discussed above. This mechanism of multiplicative coalescence, whereby large clusters grow proportionately more rapidly than small clusters, leads to the emergence of one unique giant component as discussed in more depth in Grimmett (2010); Spencer (2010). For more discussion on kernels, the gelation transition and connections to “mean-field” probability theory of random graphs see Aldous (1999); Ziff (1980). Also see Stockmayer (1943) for the kinetic theory approach to branched polymers which is equivalent to percolation on the Bethe lattice.

Note that the cluster aggregation approach, Eq. 3, assumes that each added edge connects two previously disjoint components, thus a cluster aggregation process necessarily ends at $t=(N-1)/N$ when only one component remains. If a giant component only appears at this extreme limit of $t=1$, in the limit $N\rightarrow\infty$, then this is considered lack of gelation, or trivial gelation, as there is only one cluster remaining. In contrast, alternative models of network growth beyond cluster aggregation can allow for the addition of edges that are internal to a cluster, in which case, the maximum edge density attainable would be $t=(N-1)/2$ on an undirected network. Also note that due to the equivalence of the static definition and the kinetic approach, the random graph percolation model is reversible, an issue we will return to in Sec. II.7 in the context of explosive transitions.

II.1.3 Discontinuous percolation transitions

Before the notion of Explosive Percolation (EP) was introduced in Achlioptas et al. (2009), there were many well-known studies of discontinuous percolation transitions. Several of these seminal models are summarized in Table 1. These models belong to three distinct classes. The first class is based on site or link removal, or culling, with the most prominent examples being bootstrap percolation Chalupa et al. (1979) and $k$-core percolation Seidman (1983); Bollobás (1984). An in-depth discussion of k-core percolation on complex networks can be found in Dorogovtsev et al. (2006, 2008). The second class is based on glassy dynamics Thorpe (1985). The third class is based on multiple-species models, in particular infected-susceptible type of dynamics Dodds and Watts (2004), or the 4-species models Janssen et al. (2004).

Also worth noting is the jamming transition in sphere packings as introduced in O’Hern et al. (2002). The nature of that transition was revealed by mapping the process onto analogous models related to $k$-core percolation Schwarz et al. 2006a ; Toninelli et al. (2006). These models of “jamming percolation” on low-dimensional lattices incorporate spatial correlations intended to capture glassy dynamics. They exhibit hybrid phase transitions with a discontinuous jump in an order parameter, but diverging length scales characteristic of second-order transitions Schwarz et al. 2006a ; Toninelli et al. (2006); Jeng and Schwarz (2010); Cao and Schwarz (2012).

II.2 EP from competitive percolation

Many types of graph evolution processes are now known to lead to EP. These processes delay the formation of large components and break the “multiplicative coalescence” rule of traditional models of percolation discussed in Sec. II.1.2. This allows for the many unique phenomena manifested during the emergence of large-scale connectivity, as surveyed herein. EP was first demonstrated for graph evolution processes with edge competition, which is the focus of this current Section.

II.2.1 Achlioptas Process (AP)

Random graphs are long-studied topics in combinatorics and probability theory. At a workshop entitled “Probabilistic graph theory” held at the Fields Institute in February of 2000, Dimitris Achlioptas proposed that “the power of two choices” paradigm, a recent innovation in randomized algorithms, could be used in the Erdős-Rényi graph evolution process to enhance or delay the percolation transition. This impact of two choices was first introduced in 1994 in the context of hash tables, showing that if there are two choices for the potential hash function it greatly reduces the needed storage Azar et al. (1994, 1999); Adler et al. (1998). This paradigm was shown to be particularly effective for problems such as load balancing, where giving arriving jobs a choice between the shorter of two randomly selected queues leads to exponential improvement in performance Mitzenmacher (2001). The “power of two choices” was also introduced into networking protocols based on random walks and shown to enhance performance Avin and Krishnamachari (2008). This latter model is similar to the “true” self-avoiding walk of Amit et al. (1983), where they find surprising results, such as an upper critical dimension of $d=2$.

Introducing the element of design to percolation on a lattice can also be attributed to studies of Highly Optimized Tolerance (HOT) Carlson and Doyle (1999). The intent of HOT is to use optimization, via natural selection or engineering design, to create power law distributions of cluster sizes with the clusters strategically situated to minimize damage to particular classes of random fluctuations. Specific implementations of HOT include forest fire models on a lattice. The lattice configurations achieved by the HOT process are extreme configurations that form a set of measure zero, having compact clusters that avoid a spanning cluster even at the highest densities studied. Thus, HOT models use the element of design to avoid percolation throughout the process. Occasional large failures are seen and the highly optimized states are fragile to the nature of the fluctuations.

The power of two choices for a random graph evolution process delays, but does not avoid, percolation. It can be formulated as a kinetic process as discussed in Sec. II.1.2, where edges were added uniformly at random in a discrete time manner Ben-Naim and Krapivsky (2005). But now, instead of adding one randomly chosen edge at each discrete time step, $T$, first choose two candidate edges, $\{e_{1},e_{2}\}$, and examine what would be the consequence of adding each one to the graph (an example is shown in Fig. 1(a)). Add whichever edge is more favorable by your selection criteria and discard the remaining candidate edge making it available to be a future candidate, and move on to the next discrete time step, $T+1$. Thus edges are added via a competition process that has become known as an “Achlioptas Process”. The selection criteria can vary. For instance, to enhance the onset of percolation, choose the edge that leads to the larger resulting component, and to delay the onset, choose the edge that leads to the smaller.

In 2001 Bohman and Frieze formally introduced the terminology “Achlioptas Process” (AP) and analyzed such an AP Bohman and Frieze (2001). Their goal was to use this approach to delay the formation of a giant component in an Erdős-Rényi-like graph evolution process. In order to make the model amenable to rigorous analysis, they are limited to the context of “bounded-size” rules where all components of size $K$ or greater are treated as if they were equivalent in size. For the Bohman and Frieze process (BF), $e_{1}$ is accepted if it joins two isolated nodes (and $e_{2}$ rejected), otherwise $e_{2}$ is accepted (and $e_{1}$ rejected). Thus, only components of size one (isolated nodes) are distinguished, and all components of size $K\geq 2$ are treated equivalently. A rigorous proof shows that BF delays the percolation transition when compared to ER Bohman and Frieze (2001), but this does not address the nature of the transition.

BF can be analyzed using rate equations similar to the cluster aggregation approach. Early in the evolution, well before the emergence of the giant component, all components are of size at most $O(\log N)$, and edges that are added to the graph join previously disjoint components with high probability. But in the critical regime, where components can grow sufficiently large, internal component edges become more likely to be chosen to be added to the graph. Rigorous analysis of the error introduced from violating the spanning component assumption of cluster aggregation models leads to the conjecture that all bounded-size rules lead to a continuous percolation phase transition Spencer and Wormald (2007).

II.2.2 AP with $m$-edge rules on random graphs

Achlioptas Processes with unbounded-size selection rules (where components of different sizes are treated uniquely) are not amenable to the analysis of the approaches described above Bohman and Frieze (2001); Spencer and Wormald (2007). Thus unbounded-size rules were first studied via numerical simulation of a random graph evolution process using an AP called the Product Rule (PR) Achlioptas et al. (2009).

II.2.2.1 The Product Rule and related AP’s

The Product Rule is described as follows and illustrated in Fig. 1(a). Starting from $N$ isolated nodes, two candidate edges $\{e_{1},e_{2}\}$ are chosen uniformly at random at each discrete time step, $T$. Let the vertices linked by $e_{1}$ be denoted $a$ and $b$, and the vertices linked by $e_{2}$ be denoted $c$ and $d$, and let $|C_{i}|$ denote the size of the component that contains vertex $i$. If $|C_{a}||C_{b}|<|C_{c}||C_{d}|$, then $e_{1}$ is added to the graph. Otherwise, $e_{2}$ is added. In other words, we retain the edge that minimizes the product of the two components that would be joined by that edge. Note the likelihood of choosing an edge internal to a component of size $|C_{i}|$ is proportional to $|C_{i}|^{2}/N^{2}$, which approaches zero as $N\rightarrow\infty$ in the subcritical regime when components are sub-linear in system size $N$. Thus the PR selection criteria can be used if vertices $a=b$ and/or $c=d$ without impact on the percolation transition.

Figure 1(b) shows a typical realization of a PR process, an Erdős-Rényi (ER) process Erdős and Rényi (1959, 1960); Gilbert (1959) and a Bohman-Frieze (BF) process Bohman and Frieze (2001) on a system of size $N=10^{6}$. Note that the onset of large-scale connectivity is considerably delayed for the PR process and that it emerges drastically, going from sublinear to a level approximately equal to the corresponding Erdős-Rényi and BF processes during an almost imperceptible change in edge density. Here the numerical simulations make use of the commonly used Newman-Ziff algorithm  Newman and Ziff (2001).

To quantify the abruptness of the transition, the scaling window as a function of system size $N$, denoted $\Delta_{N}(\gamma,A)$, was analyzed. This measures the number of edges required for $|C|$ to transition from being smaller than $N^{\gamma}$ to being larger than $AN$. With the specific parameters choices used in Achlioptas et al. (2009), $\gamma=A=1/2$, this measures the number of edges that need to be added for the order parameter to go from $|C|\leq\sqrt{N}$ to $|C|\geq 0.5N$. Systems up to size $N\sim 6\times 10^{7}$ were studied and the results indicated a sublinear scaling window,

[TABLE]

The associated change in edge density

[TABLE]

converging on $t_{c}\approx 0.888$ Achlioptas et al. (2009). This provides evidence of a discontinuous phase transition, but as discussed below, the Product Rule on a random graph will ultimately be shown to lead to a continuous transition in the thermodynamic limit, albeit one with an unusual universality class.

The Product Rule was then analyzed in settings beyond the Erdős-Rényi random graph, including percolation on a 2D lattice Ziff (2009) and on scale-free networks Cho et al. (2009); Radicchi and Fortunato (2009). Similar sub-linear scaling windows were observed. The scale-free networks showed additional interesting behavior. Here the degree distribution follows a power law distribution with exponent $\lambda$ and evidence suggested that there was a critical value $\lambda_{c}$ above which the system displayed a discontinuous transition, with $\lambda_{c}$ being a tricritical point Cho et al. (2009); Radicchi and Fortunato (2009). Yet, in addition, it was shown that all the models display scaling behaviors characteristic of continuous, second-order phase transitions such as power-law cluster-size distribution with an exponent close to two and scaling of susceptibility $\chi\sim|t-t_{c}|^{-\gamma}$ Cho et al. (2009); Ziff (2010); Radicchi and Fortunato (2010), but with very unusual values of scaling exponents; for a table see Radicchi and Fortunato (2010).

Many generalizations of the Product Rule and variants followed, in particular the introduction of APs with $m$-edge selection rules. Here, rather than two candidate edges, $m>2$ candidates edges are considered at each discrete time step, where $m$ is a fixed constant. For instance, in Friedman and Landsberg (2009) they introduced the “triangle rule” where the $m=3$ edges that are possible between three randomly chosen vertices are considered as candidate edges. This process was later analyzed in D’Souza and Mitzenmacher (2010). Likewise, the Sum Rule (which minimizes the sum of the components sizes that would be joined) was analyzed Cho and Kahng 2011b ; Riordan and Warnke (2012). These all show similar results of sublinear scaling windows accompanied by critical scaling behaviors. We refer the reader to the comprehensive review of these models in Bastas et al. (2014), including detailed tables of critical exponents.

II.2.2.2 The impact of a single edge

Complementary to a scaling window analysis, the impact of a single edge was studied in Refs. Nagler et al. (2011); Manna and Chatterjee (2011). This is the maximum change in the relative size of the largest component from the addition of a single edge, defined as

[TABLE]

For an $m$-edge AP, it is found that $\Delta C_{\rm max}$ decays as a power law with system size, $\Delta C_{\rm max}\sim N^{-\beta}$ Nagler et al. (2011); Manna and Chatterjee (2011). Thus, the size of this single-edge addition gap decays to zero in the thermodynamic limit. The rate of decay is typically quite small ($\beta=0.065$ for the Product Rule Manna and Chatterjee (2011)), leading to large discrete jumps in system sizes orders of magnitude larger than real-world networks. Note that the gap decaying to zero typically indicates a continuous transition, yet this also requires verifying that the growth of the order parameter displays a finite slope at the critical point as discussed formally in Sec. II.4.

A rigorous advance was made in Nagler et al. (2011) with a proof that if the largest component is allowed to grow directly in any way (i.e., by merging with another component), denoted $P_{Gr}>0$, this can lead to a continuous percolation transition. In contrast, if $P_{Gr}=0$ this necessary leads to a discontinuous percolation transition during the process. $P_{Gr}=0$ implies that growth in the size of the largest component can only happen via overtaking, when two smaller components merge together to become the new largest component. This allows for analytic calculation of the strict lower bound on the relative size of the discrete jump $\Delta C_{\rm max}\geq 1/3$ Nagler et al. (2011). The mechanism of growth by overtaking is now seen in several processes that lead to EP.

II.2.2.3 The powder keg

In Friedman and Landsberg (2009) they establish that for any percolation transition to be discontinuous there must exist a “powder keg”. This is essentially a collection of a sub-extensive number of components that together contain a total of $cN$ nodes for some constant $c>0$. Thus the number of nodes in these components diverges to infinity as $N\rightarrow\infty$. Formally the size of a powder keg, $W$, can be expressed as

[TABLE]

where $F(\tau,a)$ is the number of nodes in clusters of size greater than or equal to $a$ after the addition of the $\tau$-th edge. If $W>0$ this means there is a non-zero fraction of nodes in components ranging in size from $N^{1-\beta}$ to $N^{\alpha}$. Merging the components of such a powder keg requires only a sub-linear number of edges, leading to a discontinuous percolation transition. In fact, as shown in Friedman and Landsberg (2009), if a system is initialized with a powder keg, then even a random edge addition rule causes a discontinuous transition. For a more detailed treatment of the powder keg see Hooyberghs and Schaeybroeck (2011).

II.2.2.4 Continuous transition with unusual finite size scaling

A series of seminal works soon followed showing via numerical evidence and rate equation analysis that an AP leads, in fact, to a continuous phase transition da Costa et al. (2010); Grassberger et al. (2011); Lee et al. (2011); Tian and Shi (2012), but with many new distinct universality classes Grassberger et al. (2011); Tian and Shi (2012).

In particular in da Costa et al. (2010) they analyze a representative model of a competitive percolation process, referred to as the dCDGM model and shown in Fig. 2. Here two sets of $m$ nodes are chosen uniformly at random. The node in the smallest component of the first set is identified and then linked by an edge to the node in the smallest component of the second set. They denote the relative size of a giant component as $S$ (which is $|C|/N$ in the notation used thus far). From direct numerical simulation they show that at the critical point $S\sim(t-t_{c})^{\beta}$, indicating a continuous transition. They also develop analytic models using rate equations for the evolution in time $t$ of the size distribution $P(s,t)$ for a finite cluster of $s$ nodes to which a randomly chosen node belongs. Specifically,

[TABLE]

where $n(s)$ is the fraction of finite components of size $s$ nodes, and $\left<s\right>$ is the mean cluster size defined as the number of nodes $N$ divided by the total number of clusters that exist at time $t$. Thus, they can write the order parameter as the fraction of nodes that are not in finite components, $S(t)=1-\sum_{s}P(s,t)$.

Figure 3 shows the evolution of $P(s,t)$ for both the dCDGM model and classical random graph percolation. As shown in Fig. 3(a), although the early evolution of $P(s,t)$ for the dCDGM model (the solid black lines) suggest the build up of a powder keg, we see that at $t_{c}$ (the solid blue lines) this distribution obeys a power law, with no powder keg. Numerical solution of the rate equation for $S\approx 1-\sum_{s=1}^{10^{6}}P(s)$ leads to the value of the order parameter critical exponent $\beta=0.0555(1)\approx 1/18.$

In Grassberger et al. (2011) they also establish the continuity of APs, but using finite size scaling approaches. They denote the order parameter by $s_{\rm max}/N$ and the control parameter by $p$ (correspondingly $|C|/N$ and $t$ in this review). In particular, they focus on $P_{p,N}(s_{\rm max}/N)$ which is the distribution of the order parameter in finite systems for a given value of the control parameter. From finite size scaling theory we know that at the critical point this distribution behaves as

[TABLE]

(For details see for instance Bruce and Wilding (1992); Binder and Heermann (2010).) For a first-order (i.e., discontinuous) transition we expect $P_{p=p_{c},N}(s_{\rm max}/N)$ to be a bi-modal distribution with two distinct peaks, each one corresponding to one of the two distinct co-existing phases. We further expect that the typical distance between the peaks is the expected jump in the size of the order parameter at the critical point. In contrast, for second-order (i.e., continuous) transitions, even if the system shows the bi-modal distribution for finite $N$, we expect that as $N\rightarrow\infty$ the distance between the peaks shrinks to zero.

Shown in Fig. 4(a) is a main result from Grassberger et al. (2011), of $P_{p,N}(s_{\rm max}/N)$ for four different models of competitive percolation (the original Product Rule, the Product Rule on a 2D lattice Ziff (2009), the dCDGM model da Costa et al. (2010), and the “Adjacent edge” rule D’Souza and Mitzenmacher (2010)). The value of the control parameter $p$ is set such that the heights of both peaks are the same. They study varying $N$, extrapolating to the $N\rightarrow\infty$ limit. Although they show that all four models have continuous transitions in the thermodynamic limit, each one displays a distinct universality class, in particular values of $\beta$ and $\eta$ (see the table from Grassberger et al. (2011) reprinted in Fig 4(b)). Furthermore, all four models show double-peaked order parameter distributions, with the sharpness of the peaks increasing initially with system size. Each has different scaling laws for the width of the scaling region and for the shift of the effective $p_{c}(N)$. They conjecture that the specific non-locality of the APs studied give rise to these unusual features.

The important work of Tian and Shi (2012) also quantifies this coexistence of a strongly double peaked distribution in the histogram of the order parameter at the percolation threshold. They are able to further show that although the peaks can become more defined and separated with increasing system size, nevertheless the distance between the two peaks does shrink to zero, in a power law manner, as $N\rightarrow\infty$. They also find novel scaling exponents and universality classes. Extensive numerical studies in Lee et al. (2011) further show through a finite size scaling collapse that an AP has a well defined convergence to a continuous percolation threshold.

II.2.2.5 A rigorous proof of continuity

A rigorous mathematical proof was finally achieved in 2011 Riordan and Warnke (2011), providing analytic resolution of the $N\rightarrow\infty$ limit. In Riordan and Warnke (2011) they prove that any AP on a random graph (with a fixed number of choices) leads to a continuous percolation transition in the thermodynamic limit. In essence they show that the number of subcritical components that join together to form the emergent macroscopic-sized component is not sub-extensive in system size. There is no “powder keg”. Yet, Riordan and Warnke also showed that for an AP on a random graph, if the number of random candidate edges, $m$, is allowed to increase in any way with system size $N$, so that $m\rightarrow\infty$ as $N\rightarrow\infty$ (for instance, even as slowly as $m\sim\log(\log N))$, then this is sufficient to allow for a discontinuous percolation transition. In Waagen and D’Souza (2014), the authors show that once $m\rightarrow\infty$ as $N\rightarrow\infty$ simple selection rules based on node degree, rather than component size, can lead to EP.

In summary, rigorous analytic arguments now show that an AP on a random graph with a fixed number $m$ of candidate edges leads to a globally continuous percolation transition in the thermodynamic limit. We use the term “globally continuous” as introduced in Riordan and Warnke (2011) to mean that the order parameter is a continuous function in both the sub- and super-critical regimes. In contrast, an AP with vertex based edge selection rules (rather then edge based rules) can lack global continuity and display instead non-self averaging behaviors and discontinuous percolation transitions as discussed in Secs. II.2.4 and II.4.6. For additional details on finite-size scaling, scaling functions and critical exponents for EP and percolation more broadly, we refer the reader to a number of recent reviews Bastas et al. (2014); Araújo et al. (2014); Saberi (2015); Boccaletti et al. (2016)

II.2.2.6 AP’s to enhance percolation

As a final discussion of a standard AP on a random graph we consider enhancing (rather then delaying) the onset of percolation. For instance, consider the Product Rule but where the edge selected is the one that leads to the largest product of the sizes of the two components. This indeed leads to the earlier emergence of a giant component. But the percolation phase transition appears to have the same qualitative properties as the Erdős-Rényi model. Enhancing the onset of percolation was studied more formally in a variant that inverts the Achlioptas rule da Costa et al. (2015). The authors found a transition similar to the ordinary percolation one, though occurring in less connected systems ($p_{c}\ll 0.5$).

II.2.3 AP with $m$-edge rules on lattices

On a lattice, an AP with a fixed number $m$ of candidate edges can yield a truly discontinuous EP transition. Percolation on a lattice can be quantified by the emergence of a cluster of macroscopic size or by the emergence of a spanning cluster, the latter being a path of occupied bonds or sites that connect sites from one side of the lattice to another. Both the works by Ziff (2009) and Cho et al. (2013) study bond percolation. But, in Ziff (2009) (as discussed above) an AP is used to delay the growth of a macroscopic size cluster on a 2D square lattice. And, in Cho et al. (2013) an AP is used to delay the emergence of a spanning cluster, referred to as the spanning cluster avoidance (SCA) model, and is studied on square lattices of various dimensions $d$. The schematic of the SCA model is given in Fig. 5.

The SCA model is shown to lead to a discontinuous percolation transition for a lattice with dimension $d<d_{c}=6$ as long as the number of candidate edges for the AP is $m\geq m_{c}$, where $m_{c}$ is a constant. Here the critical value $m_{c}$ is related to $d_{BB}$, the fractal dimension of the “backbone” cluster structure which emerges close to the critical point. This is calculated analytically and measured numerically and it is shown that $m_{c}=d/(d-d_{BB})$, where $d$ is the dimension of the lattice Cho et al. (2013) . When $d=2$, $m_{c}\approx 2.554$.

An interesting distinction occurs for $m=m_{c}$ versus $m>m_{c}$. For $m=m_{c}$, the discontinuous percolation transition occurs at some intermediate $t_{c}$ during the process. In contrast, for $m>m_{c}$ the process acts globally, so when the spanning cluster emerges, it encompasses the entire system. Such global percolation also happens for an $m$-edge AP on a random graph in the limit $m\rightarrow\infty$. As first shown in in Rozenfeld et al. (2010) a giant percolating component only emerges in the final step of the process when only one component remains. This process is illustrated in Fig. 16(a). Instead of metric or geometrical confinements, the rule has unrestricted access to the entire collection of components.

II.2.4 AP with $k$-vertex rules on random graphs

The more general class of “$k$-vertex rules” (which consider a fixed number of candidate vertices rather than edges) allows for new possibilities Nagler et al. (2012); Schröder et al. (2013). Here, rather than choosing $m$ candidates edges, $k$ vertices are chosen uniformly at random and all of the $k(k-1)/2$ edges possible between them are the candidate edges. Edge selection criteria that achieve EP from $k$-vertex rules can include preference for keeping components similar in size, rather than aiming to minimize the largest component. An example of a $k$-vertex rule for $k=3$ is show in Fig. 6.

As discussed in D’Souza and Nagler (2015), this change from random-edge selection to random-vertex based edge selection ensures the ability to prohibit the direct growth of the largest component. For instance, in a $m$-edge rule, so long as there are at least two components in the system, there is a non-zero probability that all $m$ candidate edges have exactly one vertex in the largest component. It this scenario, the largest component would grow directly, thus $P_{Gr}>0$, which is compatible with a continuous transition Nagler et al. (2011). In contrast, for a $k$-vertex rule, if two or more of the $k$ vertices are in the same component, an internal component edge can be added avoiding direct growth and maintaining $P_{Gr}=0$ throughout the percolation process, which necessarily induces a discontinuous transition during the process. This process occurs in the sparse regime, where the number of edges in the graph $T\sim O(N)$, so that the probability of sampling the same edge twice decays as $1/N$.

For $k$-vertex rules the order parameter necessarily exhibits a continuous initial transition to large-scale connectivity, but a series of secondary discontinuous jumps in the order parameter are possible, with the first discontinuous jump arbitrarily close to the initial continuous transition Nagler et al. (2012); Schröder et al. (2013). This series of stochastic jumps and “devil’s staircase” transitions demonstrate an array of behaviors worthy of in-depth discussion as found in Sec. II.4.6.

II.2.5 Thermodynamic formulations of EP

In a classic work from 1972, Fortulyn and Kesten (FK) showed how to unify the framework of classic percolation on a lattice and the $q$-state Potts model Fortuin and Kasteleyn (1972). For an explicit treatment see, for instance, Grimmett (2006). In Bradde and Bianconi (2009), the authors generalize those original results of FK to extend the theory to standard percolation on complex networks. But this connection is laking for more intricate percolation problems. As discussed above, there has been intense study of the scaling exponents characterizing the critical behavior of EP processes. We refer the reader to Bastas et al. (2014) for a detailed table of critical exponents. Likewise, a Hamiltonian formulation of EP Moreira et al. (2010) is discussed in Sec. II.3, which gives insights into the underlying mechanisms of EP. In da Costa et al. (2014) the authors develop a complete scaling theory of the transition for the dCDGM model show in Fig. 2. They indicate the order parameter and the generalized susceptibility, find the full set of scaling relations and functions, including the relations between the critical exponents and the upper critical dimension.

Most recently, there have been two promising developments. First, in Bianconi (2018), the author develops a large deviation theory of percolation characterizing the response of a sparse network to rare events. The theory reveals that rare configurations of initial damage, for which the size of the giant component is suppressed, leads to a discontinuous percolation transition. A corresponding partition function is developed that allows the formulation of percolation in terms of thermodynamic quantities.

Second, the recent work of Hassan and Sabbir (2018) makes a significant advance in formulating $m$-edge Achlioptas Processes as thermal continuous phase transitions. Typically, in past work on random graphs the “susceptibility” of the random graph percolation process, denoted $\chi$, has been defined as the second moment of the component sizes as shown in Eq. (1). Instead, in Hassan and Sabbir (2018) they define susceptibility $\chi$ as a function of the ratio of successive jumps observed in the size of the largest component from addition of single edges. Although it is unexpected that such a definition would obey the fluctuation-dissipation theorem of equilibrium statistical physics, with this in place, they define thermodynamic quantities and connect explosive percolation transitions with thermal continuous phase transitions. They derive the associated scaling relationships, show that they obey the Rushbrooke inequality, and that both the Product Rule and the Sum Rule belong to the same universality class. This raises the possibility to analyze other $m$-edge AP models as thermodynamic phase transitions and may point towards the possible existence of larger universality classes in EP.

It should be noted that the results obtained in Hassan and Sabbir (2018) require that there exists a traditional scaling behavior connecting sub- and supercritical regimes. Indeed, one of the major achievements of the Wilson-Kadanoff Renormalization Group theory was to prove that such a scaling function provides a smooth connection, see for instance Wilson (1983). But, as stressed in Grassberger et al. (2011), one of the most unusual features of EP is that it seems to show completely different scalings in the sub- and supercritical regimes. So still an open question is how to resolve the results in Grassberger et al. (2011) and  Hassan and Sabbir (2018).

An important advance connecting discontinuous percolation transitions and equilibrium thermodynamics is presented in Bizhani et al. (2012). There they provide a thermodynamic formulation for a discontinuous percolation transition expressed in terms of “Hamiltonian” graphs (i.e. exponential random graph models) and, more specifically, a partition function. This allows for precise definition of thermodynamic quantities. Here chemical potentials control the number of edges, triangles and two-stars present in the graph at thermal equilibrium. They show that, for ranges of the chemical potential, the percolation transition can coincide with a first-order phase transition in the density of links, and also exhibits hysteresis loops of mixed order. Although it is a discontinuous percolation transition, it does not appear to be in the class of EP transitions which require microscopic evolution rules that delay the formation of macroscopic components. Yet the framework in Bizhani et al. (2012) may provide a promising direction for future research into the thermodynamics of explosive transitions.

II.3 EP from alternate processes

There are now several graph evolution processes known to give rise to EP transitions that do not involve direct competition between edges. These are reviewed here and reveal common mechanisms of delaying macroscopic connectivity and the common theme of the interplay of local versus global information, which is yet to be fully understood.

II.3.1 Suppressing the largest component

Several approaches achieve EP with a truly discontinuous onset of the giant component by mechanisms that directly suppress the growth of the largest component, without need for edge competition. In the “Gaussian model” of Araújo and Herrmann (2010), a regular lattice is considered for the the underlying substrate and a single edge is examined at a time (in other words $m$ = 1). If the randomly chosen edge would not increase the current size of the largest component then it is accepted. Otherwise it is rejected with a probability function that decays as a Gaussian distribution centered on the average cluster size, written as

[TABLE]

where $\langle s\rangle$ is the mean cluster size. Thus, forming components that are similar in size to the average is favored. Clear signatures of a first-order transition are observed, such as bimodal peaks for the cluster size distribution, indicating the coexistence of percolating and non-percolating regions in finite systems at $t_{c}$ as shown in Fig. 7.

Shown in Fig. 8 are examples of the clusters that form on a square lattice with $1024^{2}$ sites. In contrast to the lattice model, the equivalent random graph version of this Gaussian model exhibits a discontinuous transition at the end of the process when the system condenses into a single component Araújo and Herrmann (2010).

In Moreira et al. (2010) a Hamiltonian formalism is developed, providing a connection between equilibrium statistical mechanics and EP. They have shown that the key for obtaining a discontinuous percolation transition is that the size of the growing clusters should be kept approximately uniform. More surprisingly, they show that for random graphs, adding edges which merge together previously distinct components should dominate over adding edges internal to a component.

For more discussion of EP processes achieved by control of only the largest component, in particular by suppressing the growth of a cluster differing significantly in size from the average one, see H. J. Herrmann and N. A. M. Araújo (2011)

II.3.2 Tuning edge rejection rate (The BFW model)

A model proposed by Bohman, Frieze and Wormald (BFW) also achieves a discontinuous percolation transition by suppressing growth of the largest component and similarly considering only a single edge at a time, $m$=1 Bohman et al. (2004). Yet, the BFW model displays a broad array of phenomena worthy of discussing in its own right.

In the BFW model the system is initialized with $N$ isolated vertices, a cap on the maximum allowed component size set to $k=2$, and a budget for rejecting undesired edges where the number of edges that can be rejected increases initially with increasing $k$. In detail, a random edge is sampled from the graph, if it would not lead to an increase in $|C|$ the edge is accepted. Otherwise the edge is rejected so long as a stringent lower bound on the number of edges that must be accepted, denoted $g(k)$, is maintained. If the lower bound would be violated, instead the cap $k$ is repeatedly increased incrementally to $k+1$ until either the cap is large enough that the edge can be accepted, or $g(k)$ has decayed sufficiently that the edge can be rejected.

The formal BFW process is defined in Table 2, along with the standard notation used in the literature. So long as the fraction of accepted edges $T/u\geq g(k)$, then any candidate edge can be rejected. Here $T$ represents the number of edges accepted, and $u$ the number of discrete time steps into the process.

The lower bound $g(k)$ can be written in the most general terms as:

[TABLE]

where $\alpha$ and $\beta$ are parameters. In the original BFW model, these are not parameters, but simply constants $\alpha=1/2$ and $\beta=1/2$ Bohman et al. (2004). The system is initialized with the cap $k=2$, thus $g(k=2)=1$ and initially all edges must be accepted. The lower bound decays with increasing $k$ reaching the asymptotic limiting value of $\lim_{k\rightarrow\infty}g(k)=1/2$, meaning that asymptotically greater than or equal to one-half of all sampled edges must be accepted.

Note the tradeoff utilized by the BFW model. Early in the process when all components are small almost all edges must be accepted. But as the components grow in size, and rejection becomes more crucial for avoiding the formation of a giant component, we can reject more edges up until reaching the limiting value $g(k\rightarrow\infty)=1/2$.

The critical behavior of the BFW model was analyzed in Chen and D’Souza (2011), where the authors show that the BFW model leads to the simultaneous emergence of two co-existing giant components in a truly discontinuous percolation transition. Details are shown in Fig. 9. They also show that if edges are sampled uniformly at random from the complete graph, the co-existing giants are asymptotically stable. (Any edge that would lead to merging the two co-existing giants can always be rejected by a small sub-linear increase in the cap $k$.) If edges that are sampled are restricted to only those that join previously disjoint components, the same sub-critical behavior is observed but eventually in the supercritical regime the two components necessarily merge together causing a discontinuous jump in $|C|$.

In Chen and D’Souza (2011) the authors introduce a generalized version of the original BFW model, with the $\alpha$ parameter shown in Eq. (11). Thus the limiting fraction of accepted edges is bounded from below, not by a constant, but by a parameter $\alpha$. They show that the value of $\alpha$ determines the number of coexisting giant components that simultaneously emerge in a discontinuous transition, as shown in Fig. 10.

In Chen et al. (2012) the authors further generalize the BFW model by introducing the $\beta$ parameter shown in Eq. (11). They show that the value of $\beta$ governs the extent of direct growth allowed for the largest component. For $\beta<1$ it is always possible to reject a sampled edge that would lead to significant direct growth. Instead the growth of the size of the largest component is dominated by overtaking, when two smaller components merge together to become the new largest, giving rise to a discontinuous percolation transition. In contrast, for $\beta>1$ once the cap size $k\geq N^{1/\beta}$ there are situations when a sampled edge cannot be rejected and the largest component experiences significant direct growth, leading to a continuous percolation transition.

The BFW model displays a rich set of phenomena. Of particular interest is that $\alpha$ controls the number of co-existing giant components, and that $\beta$ controls the extent of direct growth of $|C|$ allowed. The BFW model has now been analyzed on 2D and 3D lattices, showing strong evidence of a discontinuous transition including compact clusters with fractal surfaces Schrenk et al. (2012). The supercritical properties and additional phase transitions of the BFW model on various substrates have also now been studied in for instance Chen et al. 2013c ; Chen et al. 2013b .

II.3.3 Cluster aggregation

A great number of models that display EP phenomena have been built on cluster aggregation processes, thereby ignoring the intra-cluster link structure of the underlying network. Sec. II.1.2 contains a discussion of the approach and the foundational Smoluchowski rate equation, Eq. (3), along with the definition of the collision kernel $K_{ij}$ which is proportional to the probability that components of size $i$ and $j$ merge per unit time.

The solution of the rate equation is still unknown for the great majority of collision kernels $K_{ij}$, and analytical insights are difficult. In linear polymerization, two clusters are merged by the molecules at two reactive ends, where the kernel $K_{ij}$ is a constant. When two clusters have a compact shape and merge, the kernel is given by $K_{ij}\sim(ij)^{1-1/d}$, where $d$ is the spatial dimension. A series of detailed papers based on this approach have had a deep impact on our understanding of EP systems Cho et al. (2009, 2010); Cho and Kahng 2011a ; da Costa et al. (2010); Cho and Kahng (2015); Cho et al. 2016b .

Cho et al. considered the important and more general case of a power law kernel, $K_{ij}\sim(ij)^{\omega}$, where $\omega$ is a positive exponent Cho et al. (2010). Models with $0.5<\omega\leq 1$ account for steric hindrance and intramolecular bonding Leyvraz and Tschudi (1982); Ziff et al. (1982, 1983); Leyvraz (2003). Below the critical value $\omega_{\text{c}}=0.5$ aggregation based on Eq. (3) exhibits a violation of mass conservation. In addition, for $\omega<0.5$ gelation does not take place in finite time McLeod (1962). For the modified Smoluchoswki equation with a (normalized) power law kernel

[TABLE]

total mass is conserved and gelation is guaranteed to occur in finite time Cho et al. (2010). The authors showed that for $\omega>0.5$ the transition is continuous but discontinuous for $\omega<0.5$.

As discussed earlier, controlling the growth of the largest cluster constitutes a key mechanism to delay percolation Araújo and Herrmann (2010). This motivated Cho and coworkers to study cluster aggregation where the collision rate of the largest cluster is controlled. Specifically, the exponent is taken as $\omega=\alpha$ for all reactions that do not involve the largest cluster, and otherwise $\omega=\beta$, for the general case $\alpha\neq\beta$. This model, as well as generalizations of it, show an extremely rich phenomenology of four distinct phase transition types, depending on the combination of $\alpha$ and $\beta$ (Ref.  Cho et al. 2016b ), namely continuous percolation, discontinuous percolation at the very end of the process, ultra-slow convergent discontinuous transition, and a non-self-averaging staircase behavior. We will come back to this phenomenology in Sec. II.4.6. More generally, these anomalous behaviors are expected for aggregation processes with two (or more) coalescence time scales, as artificially introduced by the two exponents Cho et al. 2016b .

In Cho and Kahng (2015) they derive the necessary conditions for the occurrence of a genuinely discontinuous phase transition for a broad class of processes that use cluster aggregation approaches. Their theory not only predicts discontinuous transitions, but can also tell the type of the transition. They showed that the key characteristic of whether or not the cluster kinetic rule is homogeneous with respect to the cluster sizes determines the phase transition type. This only requires examination of the cluster size distribution immediately before the transition. The theory is based on a mean-field two-species cluster aggregation model but may hold for a much wider range of models. As a limitation, however, the theory is expected to be non-exact for percolation processes on networks with non-random locally non-loop-free topologies, on lattices, or for other correlated or clustered underlying structures Son et al. (2011); Cai et al. (2015); Grassberger (2015); Grassberger et al. (2016).

A related model to cluster aggregation is called “agglomerative percolation” Christensen et al. (2012). Here, clusters exist on a two-dimensional lattice, but rather than diffusive motion, a discrete time process is considered. Starting from a collection of isolated clusters of size one, at each time step a random cluster is picked and bonds are added to the entire surface of the cluster in order to link it to all adjacent clusters. This is motivated in part to mimic the non-locality seen in an Achlioptas Process: If a cluster of length scale $l$ is picked, links are added simultaneously at distances $\mathcal{O}(l)$ apart. The process proceeds until the system is reduced to one cluster. Similar to multiplicative coalescence, if clusters are picked with probability proportional to their mass, this leads to a single “runaway” compact cluster. If instead all clusters are equally likely to be chosen, these leads to a continuous transition in a new universality class for the square lattice, while the transition on the triangular lattice has the same critical exponents as ordinary percolation showing intriguing violations of standard universality classes.

II.3.4 Correlated processes

The role that correlations can play in creating discontinuous percolation on a lattice has been known for some time for models of jamming, e.g., Toninelli et al. (2006); Toninelli and Biroli (2008); Jeng and Schwarz (2010). But with the growing interest in EP phenomena, in Cao and Schwarz (2012) they explore the connections between such known models and EP. In particular they analyze models that are a mixture of the correlated models and more traditional models of percolation, and they search for tricritical points separating the region of discontinuous from continuous transitions. For random graph models, they analyze a $k$-core percolation for a mixture of $k$=2-core and $k$=3-core vertices, and find a tricritical point. But for two-dimensional lattice models the behavior is not so clear, showing crossover behaviors but no tricritical point. Yet the work suggests interesting connections between EP, glassy dynamics, and jamming.

Cooperative interactions and the existence of tricritical points is also explored in Bizhani et al. (2012). In particular they analyze the tricritical point of a generalized epidemic process which maps onto a form of complex social contagion. Note several models of explosive phenomena in cooperative epidemics have also been introduced and analyzed, as discussed in-depth in Sec. V.2.

II.3.5 Hierarchical lattices

In Boettcher et al. (2012) the authors further illustrate the interplay between local and global information by examining simple percolation on a hierarchical lattice. In particular they analyze a one-dimensional lattice dressed up with a hierarchy of long-range bonds as shown in Fig. 11. It is well known that percolation on a one-dimensional lattice is trivial process with percolation only in the limit that the bond occupation probability $p=1$. Boettcher et al. recursively build up a hierarchical structure and remark that it is similar to a hyperbolic geometry where most nodes, as in a tree, are close to the periphery. They show an extensive percolation cluster arises for $p<1$ and emerges in an instantaneous jump to a finite value. Although it is not yet clear whether the transition follows the paradigm of EP phenomena, which requires microscopic dynamics that delay formation of macroscopic components, this model by Boettcher et al. provides a systematic and rigorous model to explore the interplay of local and global information that is crucial to EP but is, as of yet, still not fully explored.

II.4 Categories of EP transitions

EP transitions have been shown to display a variety of anomalous behaviors. Here we attempt to categorize the types of behaviors observed. Our motivation is to provide an umbrella that covers the main types of the critical and supercritical phenomena that have emerged from the extensive body of work on EP. It is important to note that this current classification, into only six classes, is limited and could be further expanded to incorporate more fine-grained details or as new phenomena are discovered.

As we are concerned with general universality classes, throughout this section we use the generic notation $P$ to denote the order parameter, which is the fractional size of the largest component, and $p$ to denote the control parameter.

II.4.1 Type I: Anomalous continuous

transition

We call a continuous transition that exhibits a significant macroscopic jump in size of the order parameter for any large but finite system by the term anomalous continuous, or for shorthand, we say it is of type I (see Fig. 12). A type I transition is typically characterized by the scaling of the largest gap in the order parameter as $\Delta C_{\text{max}}\sim N^{-\beta}$, for $\beta\ll 1$. This property has been used as a definition for a so-called weakly discontinuous transition D’Souza and Nagler (2015). Here we further require that for a type I transition the order parameter $P$ must have a finite slope at the critical point. This is because a gap in scaling can also be compatible with a diverging slope of $P$ (which is a combination that we call instead a type IV transition in Sec. II.4.5). As a more specific example of type I, consider some competitive rule that leads to EP and shows a weak decay of the largest gap scaling, say, $\Delta C_{\text{max}}\sim N^{-0.01}$ and no apparent divergence of slope of $P$. As discussed in D’Souza and Nagler (2015) this means that even for system sizes of the order of Avogadro’s number, $N\approx 10^{23}$, well beyond the size of real-world networks, the EP transition would still be effectively indistinguishable from a process exhibiting a genuine discontinuous transition.

Anomalous continuous EP transitions possess atypical finite-size scaling behaviors, as seen in Grassberger et al. (2011); Tian and Shi (2012), which are very different from classic continuous percolation. Once again we refer the reader to a number of recent reviews for extensive details regarding finite-size scaling, scaling functions and critical exponents for processes that lead to EP Bastas et al. (2014); Araújo et al. (2014); Saberi (2015); Boccaletti et al. (2016) and the literature therein, e.g., da Costa et al. (2014).

II.4.2 Type II: Discontinuous phase transition

A genuine discontinuous phase transition is given by a finite jump $\Delta$ of the order parameter $P$ at the percolation threshold $p_{c}<1$ as follows:

[TABLE]

where $f(p)$ is a positive function of $p$ and describes the growth in the supercritical region. Note that Eq. (13) (and the ones that will follow for the other types) hold in the thermodynamic limit. For $f(p)=(p-p_{c})^{\beta}$ and $0<\beta<1$, the discontinuous transition coincides with a second order transition (see Type IV, hybrid transitions).

Typically, however, the growth in the supercritical region, $f(p)$, is not a power law function, or if it is one, it is characterized by $\beta\geq 1$, leading to a finite slope of $P$ for $p\rightarrow p_{c}^{+}$. The required finite slope of $P$ as approached from below the critical point is indicated by the dashed line in Fig. 12.

For type II, the percolation threshold $p_{c}$ is required to be not at full link (or occupation) density, which implies that there exists an extended supercritical region. From the perspective of the kinetic approach this means that the discontinuity is required not to be at the end of the process, or in other words, gelation does occur.

Worth noting is that the majority of models of EP are not of this type II, but are either continuous, lack a supercritical region, or are of other types. Nevertheless, there are a number of models that show evidence for a type II transition on the lattice, in particular using finite size scaling Ziff (2009). However, given the challenges arising from small exponents, for the majority of models it remains open if there is a discontinuity in the thermodynamic limit, or just a weak discontinuity Ziff (2009).

For the Gaussian model Araújo and Herrmann (2010) and the spanning cluster avoidance (SCA) model, for $d<d_{c}$ and $m=m_{c}$ Cho et al. (2013), numerical evidence indicates a type II transition. As shown in Boettcher et al. (2012), ordinary percolation on a hierarchical network can be of type II, where small-world bonds are grafted onto a one-dimensional lattice.

Schröder and coworkers studied the emergence of EP from a dynamical process with link rejection on a budget specified as $b$. An efficient link rejection protocol is used to maximally delay percolation and it is assumed that rejections are costly. If the budget is small, percolation is delayed but, in contrast to other EP models, it remains in the same universality class as continuous percolation. For sufficiently large budgets, however, the percolation is maximally delayed and the transition is genuinely discontinuous Schröder et al. (2017) at $p_{c}<1$ (see Fig. 12).

The BFW model discussed in II.3.2 also implements EP by use of a “budget”. Such mechanisms break the multiplicative coalescence of classic percolation and allow for the stable coexistence of multiple giant components. As shown in detail for the BFW model, a tunable parameter $\alpha$ (the asymptotic edge rejection rate) controls the number of stable coexisting giants. Coexisting giants in the supercritical region is a novel feature of EP transitions which provides connections with modular organization and community structure in networks Rozenfeld et al. (2010); Pan et al. (2011), with more details in Sec. II.5.

II.4.3 Type III: Discontinuity at the end of the process

In many models of EP, percolation is delayed to such an extent that, in the thermodynamic limit, it takes place just at the end of the process. We assume this to coincide with full link density $p_{c}=1$, as it is on lattices and for many cluster aggregation processes (otherwise replace $p_{c}=1$ by some $p_{c}^{\text{full}}$). Thus,

[TABLE]

In the kinetic approach, full link density indicates the end of the process (as discussed in Sec. II.1.2), or equivalently that a giant component only emerges at the very end of the process and gelation does not occur.

Examples for a type III transition have been discussed in the literature Boettcher et al. (2012); D’Souza and Nagler (2015) and include 1d percolation, $m$-edge rules with infinite choice (as shown in Fig. 16(a)), cluster aggregation for power law kernels $K_{ij}\sim(ij)^{\alpha}$, for $\alpha<0.5$ (section II.3.3), and the SCA model (discussed in Sec. II.2.3) below the critical dimension $d<d_{c}$ but for sufficiently large choice $m>m_{c}$. A type III transition shows an order parameter jump at the end of the growth process (which encompasses the entire system) and lacks a supercritical region, see Fig. 12.

II.4.4 Type IV: Single-step-continuous discontinuity

In the random network model introduced by Waagen and D’Souza Waagen and D’Souza (2014), $k$ nodes are sampled at random and their respective degrees are considered. The model preferentially favors connecting nodes of small degree. If the choice parameter $k$ is an increasing function of system size, then a discontinuous transition occurs at the end of the process. Hence, the model reproduces the prediction of Riordan and Warnke regarding $k$-vertex rules for increasing $k$ Riordan and Warnke (2011). Perhaps the most remarkable particularity of the model lies in the finite size gap behavior. For increasing system size, single step gaps in the order parameter vanish in the thermodynamic limit. Yet, a genuine discontinuity emerges at the end of the process, which results from an diverging slope of $P$, see Fig. 12. Thus, the order parameter is described by Eq. (14), yet the largest gap in the relative size of the largest component $\Delta C_{\text{max}}$ decreases as a function of $N$.

A similar behavior was found for a model where multiple time scales determine irreversible cluster aggregation Cho et al. 2016b . Fig. 13 shows that the transition point is approached ultra slowly whereas the slope is an increasing function of $N$.

II.4.5 Type V: Hybrid phase transition

Random graph models displaying EP have been designed with the objective to delay percolation and typically avoid mergers of large components. The competitive percolation models discussed thus far require sampling of at least two edges, or at least three nodes in a time step. Yet, in Basatas et al. Bastas et al. (2011) and Panagiotou et al. Panagiotou et al. (2011) the authors consider “half-restricted processes” and demonstrate that sampling of only two nodes can lead to genuinely discontinuous percolation. This requires, however, one node to be sampled from a restricted set of small components, where the other node is sampled at random. Soon, a number of models based on half-restricted processes analyzed the mechanism leading to the discontinuous transition and the nature of the transition Panagiotou et al. (2013); Cho et al. 2016a ; Lee et al. (2017).

They found that the genuinely discontinuous transition coincides with a second-order transition. This combination is called a hybrid phase transition. Hybrid phase transitions exhibit both the typical critical divergence of a second-order phase transition and a genuine discontinuity,

[TABLE]

where $\Delta$ denotes the gap in the order parameter which coincides with a typical critical behavior of a power law divergence. We wish to emphasize, that a hybrid transition is a discontinuous transition with a specific class of critical behavior. Hybrid transitions require a supercritical behavior with diverging slope of the order parameter at $p_{c}$, hence $\beta<1$. Note that an exponent $\beta\geq 1$, or other supercritical behaviors with a finite slope of $P$ for $p\rightarrow p_{c}^{+}$ do not constitute a hybrid transition.

Hybrid transitions are known for models of $k$-core percolation Schwarz et al. 2006b and for cascade failure models in multiplex networks Cellai et al. (2013).

Cai et al. first reported universal mechanisms leading to hybrid transitions Cai et al. (2015). This work was followed by a study by Lee and coworkers Lee et al. (2017), where they show that a critical branching process underlies the continuous component of the transition, whereas the discontinuous component stems from supercritical branching Lee et al. (2017) with a very specific form as shown in Fig. 12. More generally, Cai and coworkers and Grassberger Cai et al. (2015); Grassberger (2015); Grassberger et al. (2016) provided compelling numerical evidence that hybrid transitions may exhibit a much larger variety of first and second order features, depending on the underlying structure or network, in particular emphasizing the role of local loops and the limitations of mean-field approximations for cascading processes Grassberger (2015).

II.4.6 Type VI: Devil’s staircases and non-self-averaging

Perhaps the most unusual feature of the models leading to EP transitions is that their behavior can be non-convergent and non-self-averaging. Likewise, continuous percolation and discontinuous percolation can coexist. These aspects were discussed briefly in Sec. II.2.4 with respect to the “Devil’s Staircase” (DS) model which is defined in Fig. 6. Here we provide the details and broader context.

Similar behaviors to the DS model are also found for the Nagler-Gutch (NG) model and the modified ER model (mER), as first reported on in Riordan and Warnke (2012). The main mechanism of both the NG and mER models is that the largest component is prevented from growing as long as the two sampled largest components are not of exactly the same size. Although this mechanism is not based on suppressing the growth of large components, it leads to anomalous percolation features. The order parameter is blurred in the supercritical regime and does not converge to a function of $p$ in the thermodynamic limit Riordan and Warnke (2012), see Fig. 12. These models show tremendous variation from one realization to another in the supercritical regime Riordan and Warnke (2012); D’Souza and Nagler (2015), see Fig. 14(a, b).

The sample-to-sample fluctuations of the NG, mER, and DS models have one “artifact” in common. For single realizations of the NG, mER, and DS model, $|C|$ necessarily jumps between the lower and upper bound envelopes Nagler et al. (2012), see Fig. 14. This rather particular behavior is a consequence of the strict impossibility for $|C|$ to grow unless the second largest cluster has exactly the same size as $|C|$, and is thus built-in to the models Riordan and Warnke (2012); Nagler et al. (2012); Schröder et al. (2013).

Note that the onset of the staircase, which is denoted as its “phase”, is determined by stochastic critical fluctuations at $p_{c}$. Yet, immediately after the largest component becomes macroscopic in size, the successive jumps in the order parameter are best described as a (deterministic) incomplete devil’s staircase. By contrast, the multi-time-scale aggregation models studied in Ref. Cho et al. 2016b can exhibit staircases that are fully stochastic. This means that not only the ensemble but also single realizations are stochastic (even for $N\rightarrow\infty$).

This anomalous behavior results from non-self-averaging, whose key property is quantified by the relative variance of the order parameter over an ensemble of realizations,

[TABLE]

where the brackets denote ensemble averaging. A system is said to be self-averaging if $\mathcal{R}_{v}(|C|)\rightarrow 0$ for $N\rightarrow\infty$. For non-self-averaging systems, however, the observable (here $|C|$) does not converge to a sharp value for large systems nor for large sample sizes. Systems that lack self-averaging therefore lack the collapse of ensemble average, minimum, and maximum as well. Non-self-averaging properties play an important role in the statistical physics of disordered systems, for instance in spin glasses Sherrington and Kirkpatrick (1975), neural networks, polymers, and population biology Derrida (1981); Mezard et al. (1987).

For second-order phase transitions it is well known that large fluctuations in $\mathcal{R}_{v}(|C|)$ are observed only in the critical window and that they collapse to a singular peak at $p_{c}$ in the thermodynamic limit.

The lack of self-averaging for the DS and the mER models is shown in Fig. 15. One observes elevated values of $\mathcal{R}_{v}(|C|)$ in the supercritical regime, which do not vanish with increasing $N$.

In the DS model, clusters that are most similar in size are preferentially merged in the process of percolation. This leads to a non-finite number of discontinuous phase transitions immediately after the first (continuous) phase transition. Single realizations for large systems are Devil’s staircases in the supercritical regime Schröder et al. (2016). Even in the thermodynamic limit, fluctuations at the critical point determine the “phase” of the staircase Schröder et al. (2013).

This stochasticity implies power law distributed crackling noise, reminiscent of the Barkhausen effect of ferromagnetism. The stochastic nature, and hence non-convergence, of the order parameter was investigated in Ref. Schröder et al. (2013) and is another example of non-self-averaging in percolation. The existence of several discrete “mini-transitions” at points that randomly fluctuate between realizations is not only true for the Barkhausen effect, but is also observed in other models of network evolution, such as the sudden emergence of “cluster cores” Foster et al. (2010), or in networks with core-periphery organization de Simón and Boguñá (2014).

II.5 EP phenomena in real-world systems

Here we review studies that have explored the paradigm of EP in real-world systems including growth of social networks, disordered materials, electric breakdown of substrates and emergence of molecular life. In many cases, the evolution of the cluster dynamics of the real system resembles those of the graph evolution models that delay formation of macroscopic components. Thus models leading to EP, and their underlying mechanisms, offer analogies and insights for real systems and allow for surprising phenomena like modular structures or non-self-averaging to arise from percolation processes. Likewise, understanding suppression mechanisms in real systems may provide future insights for the study of EP or how the paradigm might be applied to real-world systems.

II.5.1 Non-Self-Averaging and the emergence of molecular life

In 1975, Sherrington and Kirkpatrick Sherrington and Kirkpatrick (1975) introduced a mean field model for a spin glass that exhibits unusual magnetic behaviors such as non-self-averaging. Non-self-averaging has become an ubiquitous phenomenon in statistical physics of disordered systems most dominantly studied initially in spin glasses but soon the effect appeared also in neural networks, polymers, and population biology Derrida (1981); Mezard et al. (1987); Panchenko (2012). As discussed in Sec. II.4.6, non-self-averaging and phase transitions are usually defined only in the thermodynamic limit Pastur and Shcherbina (1991). In Zimmer et al. (2018), however, they argue that early molecular life is a matter of the sustained emergence of large but finite molecules and its dynamics is governed by finite size effects from the viewpoint of statistical physics. In particular, they observe anomalous behaviors of the molecule concentrations, stochastic staircases and finite-size non-self-averaging. This establishes a link between self-replication that underlies every species of living beings and an array of anomalous percolation features.

II.5.2 Finite size of real systems

In terms of the thermodynamic limit, the sizes of real-world networks, with typical ranges from millions to tens-of-billions of nodes, are quite small. The rigorous proof in Riordan and Warnke (2011) showing the lack of a powder keg implies that in the limit $N\rightarrow\infty$ the scaling window $\Delta_{N}$ is linear in system size $N$. But numerical evidence on systems up to size $N\sim 10^{7}$ indicates the window is sublinear, with $\Delta_{N}\sim N^{2/3}$ for the Product Rule Achlioptas et al. (2009). Thus, there must exist a crossover length, $N^{*}\gg 10^{7}$, where the system becomes large enough that actual realizations show convergence to the asymptotic limiting behavior.

The dCDGM model introduced in da Costa et al. (2010), illustrated in Fig. 2, enables a rough estimate of a crossover length as discussed in D’Souza and Nagler (2015). The fractional size of the largest component of the dCDGM model obeys a scaling relation $|C|/N\approx(t-t_{c})^{\beta}$ for $t$ just above $t_{c}$. For a discrete time process on a system of size $N$ the smallest increment of the control parameter possible, denoted $\Delta t$, is for one edge addition, with $\Delta t=1/N$. As the dCDGM model passes through the $t_{c}$, the next edge added causes growth

[TABLE]

We can specify a minimum jump in size of $|C|/N$ that would be observed for any finite realization of size $N$, which we denote with the parameter $0<\Delta<1$. A jump of size at least $\Delta$ is achieved for all systems of size

[TABLE]

For the dCDGM model, $1/\beta=1/0.0555\approx 18$. Setting $\Delta=0.1$ means that systems of as large as $N=10^{18}$ would show a discrete jump when passing through $t_{c}$ on the order of 10% of the system size. Thus, for the dCDMG model a crossover length requires $N^{*}>10^{18}$, far beyond the realm of most real world networks.

II.5.3 Modular network structure in the real-world

The microscopic evolution dynamics that suppress the growth of a giant component, as seen in process that lead to EP, gives rise to component structures that capture aspects of modules in a range of real-world systems.

Protein networks: An early example connecting EP and real-world processes considers an evolutionary process on the Human Protein Homology Network Rozenfeld et al. (2010). Using data on human protein homology, Rozenfeld et al. initialize each protein as disconnected and then edges between the most similar (i.e., homologous) proteins are added sequentially. (This is equivalent to an AP with global edge choice, $m=N$.) This captures the general belief that proteins evolve via duplication-mutation events from ancestral proteins, and that more similar proteins organize into network modules Medini et al. (2006). This process leads to the emergence of many large isolated components of tightly connected nodes, which eventually link together with the addition of just a few intercomponent edges, with global connectivity emerging in an explosive manner. Details are shown in Fig. 16. As the authors of Rozenfeld et al. (2010) remark, the emergent structure reflects the dense connectivity within a community and weak links between communities, suggested by Grannoveter for social systems, enabling phenomena such as the “strength of weak ties” Granovetter (1973).

Social networks: In Ref. Pan et al. (2011) the authors show that the dynamic evolution of a network following a process leading to EP can reveal underlying community structure in social networks. They analyze empirical data from two distinct real-world social networks: one is a mobile phone call network, the other is a co-authorship networks of scientists. They build up the aggregate interaction network in each case from the data and use it as the underlying substrate for a percolation process. Initially, all the edges are considered “unoccupied” and an AP is used to sequentially “occupy” edges. They show that, at $t_{c}$, the component structure reflects the underlying community structure of the network.

Global trade network: Schroder and coworkers investigate how large-scale connectivity emerges from decisions made by rational agents that individually minimize costs for satisfying their demand Schröder et al. (2018). They establish an exact mapping of the resulting nonlinear trading optimization problem to a local model that leads to EP via a large-cluster-suppression mechanism. Due to existing transportation costs, nodes satisfy their demand locally as long as it is feasible. Only at the percolation threshold are the products purchased and shipped from the globally best producer. They find that genuinely discontinuous percolation and hystereses effects can result from decision making of rational agents in the global trade network.

II.5.4 Disordered materials

Diffusion limited aggregation: The suppression of a largest component captures the evolution of several physical systems. For example, consider the seminal model of diffusion-limited cluster aggregation Witten and Sander (1981). Here clusters move via Brownian motion so that the velocity of a cluster is inversely proportional to the square root of its size and thus larger clusters move considerably more slowly. In Ref. Cho and Kahng 2011a , the authors show that diffusion-limited cluster aggregation can be mapped onto the framework of EP. They consider clusters moving on an underlying two-dimensional lattice via Brownian motion, and whenever two clusters become nearest neighbors, they merge into one larger cluster. They show that Brownian motion suppresses of the mobility of the largest clusters, impeding their growth, and leading to the discontinuous emergence of a giant cluster as a function of the number of aggregation events. They also consider a generalized Brownian motion where the velocity is inversely proportional to the mass of the cluster to a power $\eta$ and map out the tricritical point separating discontinuous from continuous emergence as a function of $\eta$.

Nanotubes: Another area where the paradigm of EP is providing new modeling opportunities is phase transitions in nanotubes. The transition from insulator to conductor is often modeled by the emergence of percolating paths in bundles of nanotubes Kim et al. (2010). But, in Ref. Kim et al. (2010), the authors show that processes leading to EP offer more realistic models as experimental studies show that the sizes of the nanotube bundles are uniform. Similar to the subcritical dynamics for EP (and unlike regular percolation), the growth of larger bundles is suppressed and the transition becomes extremely abrupt. The transition shows hysteresis, as is expected for first order transitions Kim et al. (2010).

Explosive electric breakdowns: An exciting experimental realization of explosive percolation exhibiting macroscopic electric breakdowns is a system of highly conducting particles that are adsorbed and desorbed on an insulating substrate Oliveira et al. (2014). This theoretical study links abrupt macroscopic electric breakdowns to genuine discontinuous percolation transitions.

Crackling noise: Schröder and coworkers generalized the devil’s staircase model where the merging of components with substantially different sizes are systematically suppressed and merging those whose size ratio is close to a fixed target ratio $f$ is enhanced Schröder et al. (2013). For any given target ratio $f$ a series of multiple discontinuous jumps in the supercritical regime is observed. The sizes and locations of the jumps are randomly distributed, similar to crackling noise observed in materials, such as when a sheet of paper is crumpled. This framework links non-self-averaging explosive percolation with phenomena that exhibit crackling noise, power-law fluctuations, and Barkhausen noise in ferromagnets.

Multiple-time-scale aggregation: Anomalous supercritical behaviors are expected in percolation and cluster aggregation with a separation of the reaction time scales, in particular due to mass segregation Cho et al. 2016b . Examples include systems where a force opposes diffusion in one spatial direction, as for aggregation processes affected by gravity, or in rotating (planetary ring) systems such as Saturn’s rings coalescence dynamics Brilliantov et al. (2015).

II.5.5 Explosive Immunization

The focus of this review thus far has been on delaying the onset of large-scale connectivity. Yet an AP can also be used to efficiently destroy network connectivity once it exists. In Clusella et al. (2016), they introduce “Explosive Immunization” for targeted destruction of a network on which is a disease is spreading. They use an AP to identify “superblockers”, nodes who’s removal is most efficient at destroying the network, and immunize these nodes to reduce the possibility of large epidemic outbreaks. Two different selection criteria are used depending upon whether the fraction of vaccinated nodes is sub- our supercritical, and the number of candidate edges for the AP is $m=1000$. The time complexity of the model is linear in the number of nodes $N$ (up to logarithmic corrections) allowing for efficient implementation. Testing the method on many model networks and datasets of real-world networks suggests that explosive immunization is one of the most effective strategies that exist for curtailing the spread of large-scale infection, which may prove applicable in real-world scenarios.

II.5.6 Optimal Dismantling

Dismantling of a network considers efficient node removal processes that fragmentize a given network as much as possible, subject to a fixed number of nodes or other cost functions Braunstein et al. (2016). In contrast, the problem of how to efficiently remove edges for the same objective is usually referred to as an interdiction problem (see, e.g., Ref. R. Zenklusen (2014) and references therein). Both dismantling and interdiction processes typically lead to discontinuous percolation transitions. Dismantling algorithms based on reverse EP Zhao et al. (2018) belong to the most efficient methods studied Marone and Makse (2016); Braunstein et al. (2016); Zdeborová et al. (2016).

II.5.7 Information spreading

In Zhu et al. (2018) they use the paradigm of EP to model outbreaks of information sharing in social networks, arguing that abrupt and large-scale bursts of information sharing, especially about sensitive topics, can lead to social instability. This is particularly the case if the information being spread is false. In Zhu et al. (2018), they combine models of epidemic spreading on an $L\times L$ lattice with selection rules used to determine to which neighboring node a node spreads its information. The hypothesis is that when the authenticity of the information cannot be verified, a node chooses to spread the information to the smallest group it has access to. Using the Sum Rule to implement this process in simulations, they show this approach increases the latency period when the information stays localized in small groups but at a critical stage in the evolution a global outbreak abruptly occurs, and the information is globally shared.

II.6 Other interesting EP phenomena

II.6.1 Node arrival mitigates abruptness

Two distinct studies now show that node arrival can mitigate the abruptness of an AP percolation transition, yet while still significantly delaying the onset of percolation Vijayaraghavan et al. (2013); Yi et al. (2013). In the models of percolation discussed thus far, the number of nodes $N$ is a fixed quantity and the graph evolves via edge arrival. In seminal work appearing in 2001, the impact of node arrival on the Erdős-Rényi process was analyzed Callaway et al. (2001). Starting from a few seed nodes, a new node arrives at each discrete time step and, with probability $\delta\leq 1$, an edge selected uniformly at random is added to the graph. This leads to a infinite order percolation transition similar to the Berezinskii-Kosterlitz-Thouless phase transition in condensed matter physics Callaway et al. (2001).

The random graph model of Callaway et al. (2001) is extended in Vijayaraghavan et al. (2013), so that when an edge is added to the graph (which happens with probability $\delta$), the “Adjacent Edge” rule D’Souza and Mitzenmacher (2010) is used for selecting between two candidate edges. Thus, an $m=2$-choice AP is used for edge selection, instead of choosing an edge uniformly at random. The cluster aggregation analysis in Vijayaraghavan et al. (2013) shows that this process considerably delays the onset of the percolation transition but retains the smooth, infinite order transition of the original Callaway et al. (2001) model.

In Yi et al. (2013) the authors study a similar variant of the growing random graph model of Callaway et al. (2001). Here, when it is time to add an edge to the graph, rather than adding a randomly selected edge, instead the dCDGM rule da Costa et al. (2010) for edge selection is used. Using the power-law behavior of the order parameter and analysis of the crossing of the fourth-order cumulant at the critical point, they similarly establish that the percolation transition is significantly delayed. However for this model, although the transition is smooth, it does not seem to be an infinite order transition. So the Adjacent Edge rule of Vijayaraghavan et al. (2013) and the dCDGM model Yi et al. (2013) both mitigate the abruptness of the transition, but lead to different universality classes.

II.6.2 APs with directed edges

Many important real-world networks have directed edges, where each edge points from one node to another. For instance, consider a food web where a directed edge from node $a$ to node $b$ means that species $b$ consumes species $a$. Likewise we see directed edges in gene regulation networks, scientific citation networks, the world-wide web, follower relations in on-line social networks, and many others. Yet, up until now we have only discussed networks with undirected edges.

In Squires et al. (2013) the authors explore a generalization of an AP to directed networks and study the scaling properties of the process. Many of the features of an AP on an undirected network are observed, but are less pronounced, so Squires et al. (2013) refer to AP on directed networks as “weakly explosive”. (Note in Manna and Chatterjee (2011) the authors also introduce this term “weakly explosive”, but use it to refer to the impact of a single edge discussed in Sec. II.2.2)

Of course on a directed network defining a “component” is more complicated as each node participates in an in-component (of paths to node $i$), an out-component (of paths to other nodes starting from node $i$), and potentially participate in a strongly connected component. At the percolation point for an Erdős-Rényi random graph with directed edges a giant strongly connected component, giant in-component, and giant out-component all emerge simultaneously. (See for instance Newman et al. (2001)).

In Squires et al. (2013) the authors introduce an AP (referred to as a “directed competition process”) which minimizes the product of the in-component of the node located at the tail of the edge and the out-component of the node located at the head of the edge. They analyze the critical exponents for this directed process when compared to the undirected process. In particular they find the order parameter scaling exponent $\beta\approx 0.086$ for the undirected AP but for the directed analogue the giant out-component has $\beta\approx 0.34$ and the giant strongly connected component has $\beta\approx 1.2$.

Directed networks and APs are also studied in Waagen et al. (2017). But, in addition, they consider that nodes have a rank ordering with edges preferentially directed from lower to higher-ranked individuals. This captures an essential feature of on-line social networks such as Twitter, where less popular individuals tend to follow more popular individuals (but not vice-versa). They introduce a simple percolation process (with no edge choice) on an ordered, directed network. As edges are added between disconnected nodes, it is directed monotonically with respect to the rank ordering. Large-scale connectivity occurs at very high density compared to standard percolation processes, for both the strongly connected and the weakly connected component structure, with a dominant strongly connected component appearing in a seemingly discontinuous transition. Finally, by adding a competitive percolation rule with a small bias to link users of similar rank, they show this leads to formation of two distinct components, one of high-ranked users, and one of low-ranked users, with little flow between the two components. Thus information flow is segregated into two-classes of users.

II.6.3 Degree-based rules

As discussed in Sec. II.2.2, in the seminal work Riordan and Warnke (2011) they show that for an AP on a random graph, if the number of random candidate edges, $m$, is allowed to increase in any way with system size $N$, so that $m\rightarrow\infty$ as $N\rightarrow\infty$ then this is sufficient to allow for a discontinuous percolation transition. In Waagen and D’Souza (2014) the authors show that given enough edge choice, simple local rules show EP. More specifically, they define $m(N)$ to be function describing how $m$ increases with $N$, and consider a purely local AP rule that uses only the degree of the nodes (rather than the sizes of components). They explicitly calculate the critical point, showing that $t_{c}=\left(1-\Theta(\frac{1}{m})\right)N$ and that the critical widow is of size $O\left(\frac{N}{m(N)}\right)$. Thus the scaling window is sublinear if $m$ depends on $N$, for instance $m=\log N$. Any finite realization of the process shows a stochastic jump size Waagen and D’Souza (2014) as discussed in more depth in Sec. II.4.4, exemplifying the Type IV transition in Fig. 12. They also present complementary arguments to to those in Spencer and Wormald (2007) that APs with bounded size rules lead to continuous percolation transitions.

Degree based rules are of great interest because they require only local knowledge of the node degree, whereas cluster-based rules require knowledge of a more extensive property, namely the component sizes. Recent research has focused on amplifying the effects of a local rule with the “Degree-base Product Rule” (DPR) Trevelyan et al. (2018). Here $m$ candidate edges are chosen and the one that minimizes the product of the degrees of the two nodes that would be joined is selected. Remarkably, for DPR increasing the number of choices from an initial $m=2$ does not increase the abruptness of the transition but does increasingly delay percolation. The locality of the information used constrains the evolution of the degree distribution rather than creating a powder keg of the component size distribution. Similar to node arrival Vijayaraghavan et al. (2013); Yi et al. (2013), this effectively delays the onset of percolation while mitigating the abruptness of the transition. In the limit of global information, $m=N$, DPR reduces to a model with infinite choice exemplifying the Type III transitions shown in Fig. 12.

In Hooyberghs et al. (2014) they also consider degree-based network growth, but without any competition. They consider starting from a fixed number of isolated nodes, with edges added using a degree-dependent linking probability, $p_{k}\propto k^{-\alpha}$, where $k$ is the degree of the node that would gain the edge. They can analyze the model using rate equations and show that for $\alpha>0$ nodes with low degrees are connected preferentially asymptotically approaching an EP transition. (For negative $\alpha$ incoming links are most likely attached to high-degree nodes.)

II.6.4 Explosive Ising from an AP

In Angst et al. (2012) the authors show that the concept of explosive cluster growth can also be applied successfully to kinetic spin models with cluster dynamics. They show that the use of an AP process for cluster dynamics turns the continuous Ising transition into a discontinuous transition. Explicitly they consider the two-dimensional Ising model and a modification of Swendsen-Wang dynamics. In the traditional Swendsen-Wang cluster dynamics, bonds are occupied between neighboring spins of the same orientation with probability $p$. Here, rather than occupying bonds at random, they use the Product Rule to choose between two candidate bonds. This leads to a discontinuous phase transition in magnetization that shows corresponding hysteresis behavior shown in Fig. 17. Interestingly, they demonstrate that the modification to the dynamics is expected to break detailed balance, thus the stationary state of such an explosive Ising model will no longer be an equilibrium state.

II.6.5 Discrete Scale Invariance

As discussed in Sec. II.4.6, the relative variance of EP transitions (Eq. (16)) can show exotic behaviors. It is shown in Chen et al. (2014) that peaks in the relative variance can announce the location of the critical point. In particular, the generalized BFW process (see Sec. II.3.2 and Eq. (11)) exhibits peaks in relative variance at well-defined values of link density $p_{i}$, where $i$ is an integer, which survive in the thermodynamic limit. The peaks then become infinitely sharp and their positions, $p_{i}$, follow a discrete scale invariance (DSI) pattern Sornette (1998). This pattern is described by a scale invariance that holds only for a discrete set of magnification factors. As shown in Chen et al. (2014) the critical point can be predicted through $p_{i}\rightarrow p_{c}$ from the discrete scaling relation. The non-self-averaging behaviors displayed by some EP transitions can thus provide a predictive tool for determining the location of phase transitions in these models of graph evolution. How this may be extended to real-world systems is an open problem.

II.7 General insights

II.7.1 Mechanisms and processes underlying EP

Many types of processes have now been identified that implement the delay in formation of macroscopic components necessary for EP transitions. A unifying feature of these processes is that they break the mechanism of “multiplicative coalescence” (see Sec. II.1.2) that describes standard random graph percolation.

We know from rigorous analysis that if the largest component is strictly prohibited from direct growth, meaning $P_{Gr}=0$, this necessarily leads to growth by overtaking and the discontinuous emergence of a giant component. Furthermore the minimum discontinuous jump in relative size is $\Delta C_{\rm max}\geq 1/3$ Nagler et al. (2011). The BFW process Bohman et al. (2004) comes close to achieving this mechanism. It was found that the largest component can grow directly in BFW by merging with isolated nodes, but any significant growth is from overtaking Chen et al. (2012). It is also an important observation that $k$-vertex edge choice APs can achieve $P_{Gr}=0$ but that $m$-edge choice APs cannot. Models that achieve $P_{Gr}=0$ show additional discrete jumps in $|C|/N$ in the supercritical regime after the first percolation transition Nagler et al. (2012); Schröder et al. (2013); Chen et al. 2013c ; Chen et al. 2013b .

Suppressing the growth of the largest component is a key ingredient of many processes that lead to EP. The BFW model, similar to the model of Schröder et al. (2017), achieves this by using a “budget” for edge rejection. The Gaussian model of Araújo and Herrmann (2010) suppresses direct growth in favor of creating components that are similar in size to the average component size. More generally, as shown in Moreira et al. (2010), dynamical evolution processes the keep clusters similar in size as they grow is key for obtaining a discontinuous percolation transition. The Devil’s Staircase rule Nagler et al. (2012) is an example of such a process.

The creation of a “powder keg” underlies many models that lead to EP Friedman and Landsberg (2009). Yet, finite size effects complicate its detection. For $m$-edge APs the evolution of the cluster size distribution in the subcritical regime seems to indicate the build-up of a powder keg, yet it disappears as approaching the critical point as shown in Fig. 3. The BFW process, in contrast, seems to create a true powder keg, as shown in Fig. 9.

There is also an interplay between local and global information that underlies many processes that result in EP. In general we find, from the examples reviewed here, that the more non-local the dynamical evolution rule, the easier it is to achieve EP, in that simpler rules enable EP or more abrupt transitions result.

II.7.2 The role of order parameter

In statistical mechanics the choice of the order parameter is crucial and often non-trivial Stanley (1971). In percolation, traditionally, the order parameter is the fractional size of the largest connected component of the system, giving rise to a continuous phase transition from small scale components to a unique giant component.

Many models of explosive percolation have adopted this choice of the order parameter, even if multiple macroscopic components emerge. The impact that a particular choice of the order parameter may have has been an important and controversial issue. Whenever possible, the control and order parameters should be chosen to be conjugate to each other, which is a standard practice in (canonical) statistical mechanics and thermodynamics. In percolation, however, a canonical approach is only possible for very specific models such as Hamiltonian graph models Bizhani et al. (2012).

If the order parameter is taken as the total relative size of all the macroscopic components, models of EP with an initial continuous transition but with multiple giant components, (such as the Type VI models discussed in Sec. II.4.6), may be globally continuous. This perspective neglects mesoscopic structures, yet shows that many of the observed discontinuities and non-self-averaging behaviors may be a consequence of the choice of the order parameter. In Ref. da Costa et al. (2014) the authors argued that “for continuous phase transitions the order parameter cannot be chosen in an arbitrary way, by demanding only that it is 0 in the normal phase and nonzero in the ordered phase. For these transitions the order parameter must satisfy several strict conditions that are well known in the theory of phase transitions”. In order to satisfy a number of constraints such as the hyperscaling relation, however, the order parameter needs to be based on the model details (e.g., the choice for parameter $m$) da Costa et al. (2014), even for similar models such as related $m$-edge processes.

Likewise, the choice of control parameter crucially determines which features are measured. Traditionally in percolation, the link density $p$ is taken as the control parameter. If time is used instead, for reversible percolation processes, the link density $p$ and time $t$ are equivalent. In particular, for the time-discrete case, a step forward means adding a link, whereas a step backwards means removing a link. As links are chosen entirely at random, the process of standard percolation is reversible and can be equivalently described as a reversible kinetic approach.

In EP the matter is much more intriguing because the models (continuous or discontinuous) are usually irreversible. The diffusion-limited cluster aggregation model of Ref. Cho et al. (2010) (cf. Sec. II.3.3) can exhibit a discontinuous phase transition (type III) when the order parameter is studied as a function of $p$ (the relative number of aggregation events) but not as a function of $t$, defined as the physical time of the Smoluchowski equation. Thus, choosing, from the perspective of statistical mechanics, an atypical control parameter, is sometimes delicate and can crucially alter the critical behavior and phase transition type.

Another example includes cascades in interdependent networks Buldyrev et al. (2010), which can give rise to discontinuous breakdowns. In many models, failure of a single node may cause many or even a macroscopic fraction of all nodes to fail as well. In that case, a large number of links may be removed in a single “time” step. If the largest connected component is studied as a function of time, cascading failures may show discontinuous dismantling of the network. If in each time step only a single link can be removed (as in standard percolation), cascades may show a continuous transition in the thermodynamic limit. In vivid debates it has been argued whether the link removal process occurs instantaneously with node removal, or whether the link removal processes may belong to a substantially slower time scale than the node failure time scale. In addition, Son and coworkers demonstrated that the percolation is actually continuous for fully interdependent diluted lattices. Moreover, making networks interdependent can actually smoothen the transition Son et al. (2011).

Many processes based on link culling assume instantaneous multi-link removal. Examples are numerous and include the prototypical example of discontinuous $k$-core percolation Seidman (1983); Bollobás (1984); Schwarz et al. 2006a ; Cellai et al. (2011). Generally, changing the control parameter to the relative number of removed links, thereby excluding multiple-link culling in a single step, can smoothen the transition in such a way that it becomes continuous in thermodynamic limit.

II.7.3 Reversible versus irreversible transitions

Traditionally, processes that exhibit a second order phase transition are reversible whereas first-order phase transitions generically emerge from irreversible processes. When transitions are reversible, it makes no difference whether the control parameter is tuned forwards or backwards, the order parameter shows the same continuous transition at the same critical point. Standard percolation, and models in this universality class, are of that kind.

We would like to emphasize that the occurrence of a continuous phase transition does not necessarily imply reversibility. Examples include a great majority of EP models. Typically, for explosive percolation on random networks it is impossible to construct the reversed process from the forward link-addition process, as the microscopic rule is based on cluster information which is lost during the process. For the same reason the forward process for cluster aggregation is also known to be statistically irreversible. In principle, it is possible to construct statistically reversible cluster aggregation models Kolb (1986), for example when linear molecules merge and break up independently of their size. Likewise, percolation in one dimension is an example of a reversible process that shows a discontinuous transition of type III.

Discontinuous phase transitions are usually hysteretic and often show other interesting history-dependent phenomena. Thus it is helpful to consider here an example for a reversible discontinuous percolation model in more detail. Consider the process of successively merging the two smallest clusters in the system Rozenfeld et al. (2010); Nagler et al. (2011) as shown in Fig. 16(a). For the system size $N=2^{n}$, where $n$ is an integer, the order parameter, the size of the largest component, is given by $C_{1}(T)=\frac{N}{N-T}$, where $T$ denotes the number of added links (or MC steps). This global rule generates the following sequence of the largest cluster size as 1, $\underbrace{2,\ldots,2}_{N/2\;\text{times}}$, $\underbrace{4,\ldots,4}_{N/4\;\text{times}}$, $\ldots$, with the last steps $C_{1}(N-3)=N/2$, $C_{1}(N-2)=N/2$, and finally $C_{1}=N$ at $T=N-1$. Clearly, this deterministic process can be microscopically reversed: Successively split up the largest cluster in two equal parts, which leads to the reverse sequence: $N$, $N/2$, $N/2$, $N/4$, $N/4$, $N/4$, $N/4$, $\ldots$, $\underbrace{2,\ldots,2}_{N/2\;\text{times}}$, 1. In the thermodynamic limit, this process leads to a discontinuous jump from zero to one at $p_{c}=1$ Nagler et al. (2011).

Schröder et al. discuss a special form of random fragmentation Schröder et al. (2013), where, if possible, clusters are fragmented exactly into even parts (if the cluster size is not a power of two, it is split up into clusters of size $(s+1)/2$ and $(s-1)/2$, respectively). For large system size, this process becomes the backward process of the Devil’s staircase (DS) model where only similarly sized components merge Nagler et al. (2012). Yet, as the order parameter does not converge to a function in the thermodynamic limit (Sec. II.4.6), the DS model should not be considered reversible.

III Explosive Synchronization (ES)

Explosive Synchronization (ES) describes the abrupt onset of collective synchronized behavior in a network of oscillators which arises from the interplay of the network architecture and the distribution in the natural frequencies of the oscillators. The interplay conspires to impact the threshold for the onset of synchronization and has been shown to have dynamics with hysteresis. This section introduces the standard process and variations of ES phenomenon in complex networks. In contrast to EP, which focuses on network growth models from statistical physics, ES is grounded in systems of dynamical equations. Most of the analytic studies are based on the seminal and widely studied non-linear phase oscillator model, the Kuramoto model (KM) Kuramoto (1975, 1984). This has allowed for the development of new tools and insights about the mechanisms underlying the ES process in many distinct cases. We present the detailed formulations and synthesize the classes of frameworks that lead to ES, summarizing the classes in Fig. 35 in analogy to Fig. 12 showing the classes of transitions for EP.

In addition to analytic formulations, experimental platforms Leyva et al. (2012); Kumar et al. (2015) have now implemented ES in real-world systems, allowing the confirmation of several theoretical predictions. Furthermore, in light of the developments of ES, many abrupt transitions and hypersensitive responses to perturbations observed in neuroscience have been recently revisited from the perspective of ES Kim et al. (2016, 2017); Wang et al. (2017, 2016); Lee et al. (2018), and the results suggest ES can be a useful paradigm for understanding neural phenomena from epilepsy to the conscious-unconscious transition when awaking from anesthesia.

Note that discontinuous synchronization transitions are known to show up also in mean field models but arising from mechanisms different to those discussed here. For instance, in the mean field Kuramoto model, Tanaka et al. 1997a ; Tanaka et al. 1997b showed that adding inertial terms can lead to a discontinuous transition and Pazo showed that a discontinuous transition emerges when the frequency distribution becomes sufficiently broad Pazó (2005). However, explosive network transitions are rooted in mechanisms that involve both the local connectivity among nodes and the dynamical process. In particular, from the diverse mechanisms leading to ES that are summarized later in Fig. 35, it is easy to note that most of the corresponding ES frameworks have no counterpart in mean-field models since their application demands the existence of a networked structure connecting dynamical units.

III.1 Introduction to the Kuramoto model

Before presenting the recent advances on ES phenomena starting in Sec. III.3, let us briefly review here the basics of the Kuramoto model (KM) which is the dynamical model analyzed in most of the ES studies. The KM implemented on networks of coupled units is reviewed in Sec. III.2. Collective synchronization arises when a large number of interacting units oscillate in phase together uniformly regardless of the differences in their intrinsic dynamics. This phenomenon has attracted the attention of scientists for quite some time due to its beauty and its widespread presence in many diverse fields of science. This ubiquity ranges from networks of pacemaker cells in the heart and flashing fireflies to many examples in physics and engineering such as arrays of lasers and superconducting Josephson junctions Strogatz (2000); Strogatz and Goldenfeld (2004); Acebrón et al. (2005). The pioneering work by Winfree (1967) on collective synchronization spurred the field and called for mathematical approaches to tackle the problem. One of these approaches considers that the system is made up of a huge population of coupled, nearly identical, interacting limit-cycle oscillators, where each oscillator exerts a phase dependent influence on the others and changes its rhythm according to a coupling function Strogatz (2000); Acebrón et al. (2005).

Even if the previous simplifications seem to be very crude, the phenomenology of the synchronization problem can be obtained. When the natural (intrinsic) frequencies of the oscillators are too diverse compared to the strength of the coupling, they are unable to synchronize and the system behaves incoherently, showing a range of phases whose average is nearly zero. However, if the coupling is strong enough, the behavior of all the oscillators aligns and they beat in synchrony. The transition from incoherent to coherent behavior takes place smoothly at a certain threshold where synchrony starts to gradually emerge. At the transition point some set of elements lock their relative phase forming a cluster of synchronized nodes marking the onset of synchronization. Beyond this value, the population of oscillators is split into a partially synchronized population made up of oscillators locked in phase and a group of other oscillators whose natural frequencies are too different to allow them to join the coherent cluster. Finally, after further increasing the coupling, more and more elements get entrained around the mean phase of the collective rhythm generated by the whole population and the system settles into the completely synchronized state, see Fig. 18.

Kuramoto Kuramoto (1975, 1984) worked out a mathematically tractable model to describe this phenomenology, which we denote the Kuramoto Model (KM). He noticed that the most suitable case for analytical treatment should be the mean field approach, which corresponds to equally weighted, all-to-all coupling, taken to be purely sinusoidal for simplicity. In particular, Kuramoto proposed the governing equations for each of the $N$ oscillators $i$ in the system:

[TABLE]

where the factor $1/N$ is included to avoid divergences in the thermodynamic limit ($N\rightarrow\infty$), $\omega_{i}$ stands for the natural frequency of oscillator $i$, and $K$ is the coupling constant. The frequencies $\omega_{i}$ are distributed according to some function $g(\omega)$, that is usually assumed to be unimodal and symmetric about its mean frequency $\Omega$. The most commonly used frequency distributions are the Lorentzian (Cauchy) ($g(\omega)=\frac{1}{\pi}\frac{\gamma}{(\omega-\Omega)^{2}+\gamma^{2}}$) and Gaussian ($g(\omega)=\frac{1}{\sqrt{2\pi\sigma}}{\mbox{exp}}\left(-\frac{(\omega-\Omega)^{2}}{2\sigma^{2}}\right)$) distributions. Due to the rotational symmetry in the model, we can use a rotating frame and redefine $\omega_{i}\rightarrow\omega_{i}+\Omega$ for all $i$ and set $\Omega=0$, so that the $\omega_{i}$’s denote deviations from the mean frequency.

The whole population’s collective dynamics is measured by the macroscopic complex order parameter,

[TABLE]

where the modulus $0\leq r(t)\leq 1$ measures the phase coherence of the population and $\phi(t)$ is the average phase. The values $r\simeq 1$ and $r\simeq 0$ (where $\simeq$ stands for fluctuations of size $O(N^{-1/2})$) describe the respective limits in which all oscillators are phase locked and in which they all move incoherently. Multiplying both parts of Eq. (19) by $\text{e}^{\text{i}\theta_{i}}$ and substituting the imaginary part into Eq. (19) yields

[TABLE]

Equation (21) states that each oscillator interacts with all the others only through the mean field quantities $r$ and $\phi$. The second term provides a positive feedback loop to the system’s collective rhythm. Namely, as $r$ increases the population becomes more coherent, so the coupling between the oscillators is further strengthened and more of them can in-turn be recruited to take part in the coherent set. Moreover, Eq. (21) allows for the calculation of the critical coupling parameter $K_{c}$ and the characterization of the order parameter $\lim_{t\rightarrow\infty}r_{t}(K)=r(K).$ Looking for steady-state solutions, one assumes that $r(t)$ is constant and $\phi(t)$ rotates uniformly at frequency $\Omega$. Without loss of generality we can set in Eq. (21) $\phi=0$ and go into the rotating frame Kuramoto (1975, 1984); Acebrón et al. (2005). This reveals two different long-term behaviors when the coupling is larger than the critical value, $K_{c}$. On the one hand, when $|\omega_{i}|\leq Kr$, a group of oscillators are phase-locked at frequency $\Omega$ in the original frame according to the equation $\omega_{i}=Kr\sin{(\theta_{i})}$. On the other hand, the rest of oscillators for which $|\omega_{i}|>Kr$ is verified, are drifting around the circle, sometimes accelerating and sometimes rotating at lower frequencies. Demanding that some conditions hold for the stationary distribution of drifting oscillators with frequency $\omega_{i}$ and phases $\theta_{i}$ Acebrón et al. (2005), a self-consistent equation for $r$ can be derived as

[TABLE]

which admits a non-trivial solution corresponding to the existence of partially synchronized oscillators given by

[TABLE]

Equation  (22) corresponds to the mean field approach for the critical coupling at the onset of synchronization. Moreover, near the onset, the order parameter, $r$, obeys the usual square-root scaling law for mean field models: $r\sim(K-K_{c})^{\beta}$ with $\beta=1/2$. Numerical simulations of the model have proven that the mathematical approach developed by Kuramoto matches the numerics. However, in Coletta et al. (2017) the authors prove that the coupling range $\delta K$ over which this scaling is valid shrinks like $\delta K\sim N^{-\alpha}$ with $\alpha\approx 1.5$ as $N\rightarrow\infty$. Away from this interval, the order parameter exhibits the scaling behavior $r-r_{c}\sim(K-K_{c})^{2/3}$ predicted by Pazó (2005). In this latter work, a first-order phase transition was revealed from imposing a uniform frequency distribution on the Kuramoto model. The transition is second-order for other distributions of frequencies that are not strictly uniform (Lorentzian Gaussian, parabolic, etc.).

The KM approach to synchronization was a great breakthrough for the understanding of the emergence of synchronization in large populations of oscillators. Even in the simplest case of a mean field interaction, there remain unsolved problems that have resisted any analytical attempt. This is the case, e.g., for finite populations of oscillators with respect to some questions regarding global stability results Bonilla et al. (1992); Acebrón et al. (2005). In what follows, however, we focus on another aspect of the model’s assumptions, namely that the connection topology of real systems Newman 2003b ; Boccaletti et al. (2006) usually do not show a homogeneous pattern of interconnections as required by the mean field approach, but instead posses complex patterns encoded as a network, see reviews Arenas et al. (2008); Dörfler and Bullo (2014). Note, the KM is still a subject of relevance in understanding biological systems, for example, very recently Hannay et al. (2018) introduced a macroscopic reduction for networks of coupled oscillators motivated by a structure found in experimental measurements of circadian protein expression and in several mathematical models for coupled biological oscillators. They provide the first evaluation of the suitability of the Ott-Antonsen Ott and Antonsen (2008) reduction procedure for extracting macroscopic models of real biological networks.

III.2 The Kuramoto model in networks

The KM was also one of the first dynamical processes to be analyzed on a complex network topology. However, to do so, it was necessary first to reformulate Eq. (19) to include the network connectivity between the oscillators:

[TABLE]

where $\lambda_{ij}$ is the coupling strength between pairs of connected oscillators and $A_{ij}$ are the elements of the adjacency matrix encoding the connections between the nodes of the network. Note that the original KM is recovered by letting $A_{ij}=1,\forall i\neq j$ (all-to-all) and $\lambda_{ij}=K/N,\forall i,j$. In what follows, we will also consider $\lambda_{ij}=\lambda$, so that the coupling strength between any pair of connected nodes is identical.

Studies of synchronization on complex network topologies, where each node is considered to be a Kuramoto oscillator, were first reported for Watts-Strogatz (WS) networks Watts (1999); Hong et al. (2002) and Barabási-Albert (BA) scale-free (SF) graphs Moreno and Pacheco (2004); Moreno et al. (2004), see Rodrigues et al. (2016) for a review. These works are mainly numerical explorations of the onset of synchronization, with the main goal being identifying the critical coupling above which sets of nodes beating coherently first appear. In Hong et al. (2002), the authors considered oscillators with intrinsic frequencies distributed according to a Gaussian distribution with unit variance that are arranged in a WS network. They explored how the order parameter, Eq. (20), changes upon addition of long-range links for varying rewiring probability, $p$. Moreover, they assumed a normalized coupling strength $\lambda_{ij}=\lambda/\langle k\rangle$, where $\langle k\rangle$ is the average degree of the graph. Numerical integration of the equations of motion (23) under variation of $p$ shows that collective synchronization emerges even for very small values of the rewiring probability.

The results confirm that networks obtained from a regular ring with just a small fraction of rewired links ($p\gtrsim 0$) can be synchronized with a finite $\lambda$. Moreover, in contrast with the arguments provided in Hong et al. (2002), we notice that their results had been obtained for a fixed average degree and thus the KM’s critical coupling cannot be recovered by simply taking $p\rightarrow 1$ (which produces a random ER graph with a fixed minimum connectivity). This limit is recovered by letting $\langle k\rangle$ increase. More specifically, numerical simulations of the same model in Watts (1999) showed that the Kuramoto limit is approached as the average connectivity grows.

The existence of a critical point for the KM on SF networks came as a surprise. This is one of the few cases in which a dynamical process shows a critical behavior when the substrate is described by a power-law connectivity distribution with an exponent $\gamma\leq 3$ Newman 2003b ; Boccaletti et al. (2006); Dorogovtsev et al. (2008). In principle it could be a finite size effect, but results from numerical simulations show that this is not the case. To determine the exact value of $\lambda_{c}$, one can make use of standard finite-size scaling analysis. At least two complementary strategies have been reported. The first one allows bounding the critical point and is computationally more expensive. Consider a network of size $N$, for which no synchronization is attained below $\lambda_{c}$, where $r(t)$ decays to a small residual value of size $O(1/\sqrt{N})$. Then, the critical point may be found by examining the $N$-dependence of $r(\lambda,N)$. In the sub-critical regime ($\lambda<\lambda_{c}$), the stationary value of $r$ falls off as $N^{-1/2}$, while for $\lambda>\lambda_{c}$, the order parameter reaches a stationary value as $N\rightarrow\infty$ (though still with $O(1/\sqrt{N})$ fluctuations). Therefore, plots of $r$ versus $N$ allow us to locate the critical point $\lambda_{c}$.

Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is worth briefly summarizing other numerical results that have explored how the critical coupling depends on other topological features of the underlying SF graph. Recent results have shed light on the influence of the topology of the local interactions on the route to and the onset of synchronization. In particular, the authors in McGraw and Menzinger (2005, 2007); Gómez-Gardeñes and Moreno (2007) explored the Kuramoto dynamics on networks in which the degree distribution is kept fixed, while the clustering coefficient ($C$) and the average path length ($\ell$) of the graph change. The results suggest that the onset of synchronization is mainly determined by $C$, namely, networks with a high clustering coefficient promote synchronization at lower values of the coupling strength. On the other hand, when the coupling is increased beyond the critical point, the effect of $\ell$ dominates over $C$ and the phase diagram is smoothed out (a sort of stretching), delaying the appearance of the fully synchronized state as the average shortest path length increases.

In a series of papers Ichinomiya (2004, 2005); Restrepo et al. 2005a ; Restrepo et al. 2005b ; Lee (2005); Restrepo et al. (2006), the onset of synchronization in large networks of coupled oscillators has been analyzed from a theoretical point of view under certain (sometimes strong) assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are available up to date. Moreover, the analytical approaches predict that for uncorrelated SF networks with an exponent $\gamma\leq 3$, the critical coupling vanishes as $N\rightarrow\infty$, in contrast to numerical simulations on BA model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity strongly hampers analytical solutions of the model.

The analytical insights discussed so far can be formulated in terms of a mean field approximation for complex networks Ichinomiya (2004); Lee (2005); Restrepo et al. 2005a ; Ichinomiya (2005). This approach (valid for large enough $\langle k\rangle$) considers that every oscillator is influenced by the local field created in its neighborhood. To this aim, a local order parameter defined as:

[TABLE]

is introduced, where $\langle\cdot\cdot\cdot\rangle_{t}$ stands for a time average. Considering that the degree of a node $i$ is large enough, and that the network is close to the onset of synchronization, $\lambda\sim\lambda_{c}$, it is possible to write the local order parameter $r_{i}$ as:

[TABLE]

The definition of $r_{i}$ allows us to define a new global order parameter $r$ to measure the global coherence. Moreover, assuming that $r_{i}$ is proportional to the degree of the corresponding node $k_{i}$, i.e., $r_{i}\sim k_{i}$, one can write $r$ as:

[TABLE]

After summing over $i$ and substituting $r_{i}=rk_{i}$ in Eq. (25) we obtain:

[TABLE]

The above relation, Eq. (27), was independently derived in Ichinomiya (2004), who first studied analytically the problem of synchronization in complex networks, though a different approach. Taking the continuum limit and considering the limit $r\rightarrow 0^{+}$ ($\lambda\rightarrow\lambda_{c}^{+}$), Eq. (27) becomes

[TABLE]

which leads to the condition for the onset of synchronization:

[TABLE]

The mean field result, Eq. (29), gives us a surprising result that the critical coupling $\lambda_{c}$ in complex networks is nothing else but the one corresponding to the all-to-all topology $K_{c}$ re-scaled by the ratio between the first two moments of the degree distribution, regardless of the many possible differences between the patterns of interconnections of the underlying network. Precisely, it states that the critical coupling strongly depends on the topology of the underlying graph. In particular, $\lambda_{c}\rightarrow 0$ when the second moment of the distribution $\langle k^{2}\rangle$ diverges, which is the case for SF networks with $\gamma\leq 3$. Note, the contrast with the result for $\gamma>3$, when the coupling strength does not vanish in the thermodynamic limit. On the other hand, if the mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing $\gamma$, the onset of synchronization occurs at smaller values of $\lambda_{c}$. Note that Eq. (29) is mathematically correct only for locally tree-like networks, however as pointed out in Ichinomiya (2005), the approximation works well even for clustered networks, depending on the probabilities of the different link forming triangles. This surprising effect was scrutinized later to unveil the unreasonable effectiveness of tree-based theory for networks with clustering Melnik et al. (2011) where the authors proved that such a theory works well as long as the mean intervertex distance l is close to the expected value of l in a random network with negligible clustering and the same degree-degree correlations.

Interestingly enough, the dependence derived in Eq. (29) has the same functional form for the critical points of other dynamical processes such as percolation and epidemic spreading processes Newman 2003b ; Boccaletti et al. (2006); Dorogovtsev et al. (2008). While this result is still under numerical scrutiny, it would imply that the critical properties of many dynamical processes on complex networks are essentially determined by the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this last claim will be of extreme importance in physics, probably changing many preconceived ideas about the nature of dynamical phenomena.

Within the mean field theory, it is also possible to obtain the behavior of the order parameter $r$ near the transition to synchronization. Another theoretical study in Oh et al. (2007) is worth mentioning here. They have extended the mean field approach to the case in which the coupling is asymmetric and depends on the degree. In particular, they studied a system of oscillators arranged in a complex topology whose dynamics is given by

[TABLE]

Here, $\eta=1$ corresponds to the symmetric, non-degree dependent, case. Extending the mean field formalism to the cases $\eta\neq 1$, they investigated the nature of the synchronization transition as a function of the magnitude and sign of the parameter $\eta$. By exploring the whole parameter space $(\eta,\lambda)$, they found that for $\eta=0$ and SF networks with $2<\gamma<3$, a finite critical coupling $\lambda_{c}$ is recovered in sharp contrast to the non-weighted coupling case for which we already know that $\lambda_{c}=0$. This result seems phenomenologically meaningful, since setting $\eta=0$ implies that the coupling in Eq. (23) is $\lambda_{ij}=\lambda/k_{i}$, which, as discussed earlier Motter et al. (2005), may have the effect of partially destroying the heterogeneity inherent to the underlying graph by normalizing all the contributions $\sum_{j=1}^{N}A_{ij}\sin(\theta_{j}-\theta_{i})$ by $k_{i}=\sum_{j=1}^{N}A_{ij}$.

In a more complex scenario Lopes et al. (2016) the authors addressed the impact of random fields on the emergence of synchrony in the Kuramoto model on complete (fully-connected) graphs and uncorrelated random complex networks. Considering random fields with uniformly distributed directions and homogeneous and heterogeneous (Gaussian) field magnitude distribution, they prove that the critical behavior of the order parameter presents a tricritical point above which a second-order phase transition gives place to a first-order phase transition when the network is either fully connected or scale-free with the degree exponent $\gamma>5$. The tricritical point is at the intersection of three regions: (I) asynchronous state, (II) region of hysteresis, and (III) partially synchronous state. Below the tricritical point there is a second-order phase transition, whereas above it there is a first-order phase transition, see Fig. 19.

Moreover, for scale-free networks with exponents $2<\gamma<5$, the phase transition is of second-order at any field magnitude, except for degree distributions with $\gamma=3$ when the transition is of infinite order at $\lambda_{c}=0$ independent of the random fields.

III.3 Explosive Synchronization in networks

We now introduce the basics of ES phenomena as a genuine effect of the interplay between the structure of the interaction networks and oscillators’ dynamics. As already mentioned, discontinuous transitions in the usual (all-to-all) Kuramoto model were found in Pazó (2005) by using broad and homogeneous frequency distributions $g(\omega)$. However, the abrupt transitions observed in the synchronization diagram $r(\lambda)$ do not lead to any bi-stable region in which, for a range of coupling strengths, the synchronized and the incoherent states coexist. In addition, in degree heterogeneous networks the use of this kind of frequency distributions does not alter the nature of the transition. On the contrary, the use of homogeneous frequency distributions in scale-free networks Moreno and Pacheco (2004); Gómez-Gardeñes et al. 2007b corroborates through finite size scaling the second-order nature of the transition, with $\beta\sim 1/2$. This finding was obtained together with the expectation that the onset of synchronization is enhanced as degree heterogeneity increases (controlled by the second moment of the degree distribution $\langle k^{2}\rangle$) as reported in Eq. (29).

For the great majority of abrupt ES phenomena there are no rigorous proofs for the occurrence of a genuine discontinuity. Yet, there is often compelling numerical evidence for it. In other words, the phase transition categorization of ES, in analogy to that sketched-out in Fig. 12 for EP, is unknown. As a result, we will call explosive synchronization transitions which may be genuinely discontinuous by the term abrupt, to avoid unnecessary overstatements.

III.4 The correlated degree-frequency (CDF) framework

As introduced above, the ability of scale-free networks to synchronize relies on the role of hubs in centralizing the growth of a single synchronized connected component Gómez-Gardeñes et al. 2007a . Exploiting the idea that by delaying the onset of a collective behavior, such as synchronization, it is possible to change the nature of the corresponding transition, Gómez-Gardeñes et al. Gómez-Gardeñes et al. 2011a proposed a modified Kuramoto model in which the frequency of each oscillator is correlated with its degree, $\omega_{i}=f(k_{i})$. This degree-frequency correlation introduces an interplay between structural and dynamical attributes of the nodes. In particular, they chose to investigate the simplest linear correlation, $\omega_{i}=k_{i}$, thus ofsetting the topological ability of hubs to synchronization (provided by their large connectivity) by assigning them a very disparate natural frequency. Under this framework the Kuramoto evolution equations become

[TABLE]

In addition, the perfect correlation between degrees and natural frequencies automatically sets: $g(\omega)=P(k)$ and, consequently, $\langle k\rangle=\Omega$.

The numerical study of this model was done in homogeneous and scale-free networks by monitoring the synchronization transition $r(\lambda)$ though continuation, i.e. by increasing (forward) and decreasing (backward) the coupling $\lambda$, in steps of length $\delta\lambda\ll 1$. This way, the initial condition for each $\lambda$ value is exactly the final configuration, $\{\theta_{i}\}$, obtained for the previous value of the coupling strength. In Fig. 20(a) the forward and backward branches of the synchronization diagram are shown for homogeneous ER networks. In this case the transition is smooth and both branches are equal thus indicating that the transition is, as usual, of second order. In addition, the authors computed the evolution of the effective frequency of each node defined as

[TABLE]

with $T\gg 1$. The evolution with $\lambda$ of these effective frequencies is shown in Fig. 20(b) for the ER network. From this panel we observe that some of the first oscillators to converge to the average frequency ($\Omega=\langle k\rangle=6$ in this network) are those have the largest degree. This observation is confirmed by computing (see solid lines in Fig. 20(b)) the average effective frequency of nodes having the same degree (and hence the same natural frequency) as

[TABLE]

where $N_{k}=NP(k)$ is the number of nodes with degree $k$ in the network.

The above scenario dramatically changes when considering SF networks (see Figs. 20(c) and 20(d)). In this case the synchronization diagrams show an abrupt phase transitions and, moreover, the forward and backward branches do not coincide, pointing out hysteresis phenomenon. The abrupt change of $r(\lambda)$ becomes evident by looking at the evolution (by increasing $\lambda$) of the effective frequencies as defined in Eqs. (32) and (33). At variance with the ER network, the effective frequencies of all the nodes remain roughly locked to their natural ones ($\omega_{i}=k_{i}$) until the critical (forward) coupling, $\lambda\sim 1.42$, for which all the effective frequencies collapse into the same value $\Omega=\langle k\rangle=6$.

III.4.1 Analytical characterization

The numerical finding of an explosive synchronization transition in scale-free networks motivated the theoretical analysis of the Kuramoto model with degree-frequency correlations as proposed in Gómez-Gardeñes et al. 2011a . In particular, in Peron and Rodrigues 2012a the expression of the critical coupling, $\lambda_{c}$ was analyzed making use of the framework introduced in Ichinomiya (2004) for the analysis of the critical properties of the Kuramoto model in complex networks. This way, by considering the density of nodes with degree $k$ and phase $\theta$ at time $t$, $\rho(k;\theta,t)$ and taking the continuum limit for the degree variable, Eq. (31) can be rewritten as

[TABLE]

where $P(k^{\prime}|k)$ is the probability that, given a node of degree $k$, it has a neighbor of degree $k^{\prime}$. For uncorrelated networks this is simply: $P(k^{\prime}|k)=k^{\prime}P(k^{\prime})/\langle k\rangle$. After some derivations in the same line of Ichinomiya (2004), the Kuramoto order parameter, $r$, is found to satisfy the following equation

[TABLE]

The non-trivial solution $r\neq 0$ close to the synchronization onset ($r\rightarrow 0^{+}$) yields the critical coupling

[TABLE]

Note that this expression (for the threshold for the degree correlated framework) is different from the generic one for the synchronization threshold of the Kuramoto model in a network shown in Eq. (29). The validity of Eq. 36 was numerically tested for different synthetics networks.

In a later work Pinto and Saa (2015), the analytical derivation of the critical coupling was generalized to the case in which degree-frequency correlation only holds for oscillators with $k\geq k^{*}$, while for nodes with $k<k^{*}$ the frequency is randomly assigned following a distribution $g(\omega)$ as in the usual Kuramoto model. Following again the approach introduced in Ichinomiya (2004) the expression for the critical coupling becomes

[TABLE]

where

[TABLE]

It is easy to check that when $k^{*}=k_{min}$, so that the degree-frequency correlation holds for all of the oscillators, then $\alpha=0$ and $\beta=\langle k\rangle$, recovering the expression for the critical coupling shown in Eq. (36).

The dependence of the explosive transition on the degree heterogeneity of the underlying network was analyzed in detail in Countinho et al. (2013). By using an annealed network approach, i.e. working with an effective adjacency matrix with elements $A_{ij}=k_{i}k_{j}/(2\langle k\rangle)$, the authors showed that the nature of the phase transition in scale-free networks depends on the exponent $\gamma$ of the degree distribution, $P(k)\sim k^{-\gamma}$. In particular, under the assumption of the annealed approach, scale-free networks with $2<\gamma<3$ display a first order phase transition while for $\gamma>3$ the transition becomes second order. Interestingly, in the case $\gamma=3$, the CDF Kuramoto model undergoes an hybrid phase transition, i.e. the synchronization order parameter $r$ shows an abrupt change yet displays the critical properties typical of second order transitions. The authors show that for $\gamma=3$ the critical behavior close to the synchronization onset is: $r-r_{c}\propto(\lambda-\lambda_{c})^{2/3}$. Figure 21 shows the phase diagram in the $(\lambda,\gamma)$-plane. This diagram shows, for $\gamma<3$, a phase (III) where the fully synchronized state and the unsynchronized solution coexist and it is this which leads to the hysteresis phenomenon. Remarkably, the region associated with this phase increases as the exponent $\gamma$ of the degree distribution decreases, indicating that the hysteresis cycle becomes larger in this regime.

III.4.2 The star network

A common way to explore collective phenomena found in scale-free networks is to analyze star-like networks. These networks capture in a minimal way the degree heterogeneity found in scale-free networks while allowing analytical calculations without strong assumptions. In the original article showing the existence of ES in the CDF Kuramoto model Gómez-Gardeñes et al. 2011a the authors analyzed this toy graph composed of a central hub and $K$ peripheral nodes. Under the CDF scheme the frequency of the hub becomes $\omega_{h}=K$ whereas the leaves (i.e., nodes of degree one) have $\omega=1$. This way, the Kuramoto evolution becomes

[TABLE]

and the average frequency is $\Omega=2K/(K+1)$.

The numerical study of the star network is shown in Fig. 22. The top panel shows the backward and forward synchronization diagrams for a star network with $K=10$ leaves and illustrates the hysteresis cycle that is bounded by the two coupling values $\lambda_{c}^{b}$ and $\lambda_{c}^{f}$. In Gómez-Gardeñes et al. 2011a the authors derived the exact value of the critical coupling $\lambda_{c}^{b}$ by considering the backward branch and inspecting the necessary conditions for the existence of the locked solution as an attractor of the dynamics. To this aim, one can set the phase of the hub as the reference point by defining $\phi_{i}=\theta_{h}-\theta_{i}$ so that the set of Eqs. (40) and (41) transforms into

[TABLE]

In the locked regime $\dot{\phi}_{i}=0$ and $\phi_{i}=\phi$ $\forall i$ so that the condition for the synchronized solution to exist is

[TABLE]

The minimum $\lambda$ value that satisfies the former equation is

[TABLE]

Remarkably, the value $\lambda_{c}^{b}\rightarrow 1$ as $K$ increases, thus being almost independent of the number of leaves surrounding the hub. This result was also obtained in Zou et al. (2014) by analyzing the local attractiveness of the synchronized solution, allowing them to derive as well the value of the synchronization parameter $r$ at $\lambda_{c}^{b}$, $r_{c}^{b}$ (see Fig. 22), i.e. when the locking manifold disappears,

[TABLE]

An estimation of $r_{c}^{b}$ was also derived in Gómez-Gardeñes et al. 2011a yielding a slightly different result: $r_{c}^{b}=K/(K+1)$, which coincides with Eq. (45) when $K\gg 1$.

The calculation of the upper bound, $\lambda_{c}^{f}$, of the hysteresis cycle is more complicated and was tackled in Zou et al. (2014). Making use of the theory developed in Pereira et al. (2013), the authors estimate the value $\lambda_{c}^{f}$ by identifying the value of $\lambda$ for which the basin of attraction of the locked manifold passes from being local (and thus coexisting with the unlocked solution) to global. Their derivation yield a value of: $\lambda_{c}^{f}\simeq 0.6989(K-1)/\sqrt{K})$. This value was later calculated in Xu et al. (2015) yielding

[TABLE]

that for $K\gg 1$ coincides with the one derived in Zou et al. (2014). These two results point out a scaling of the forward critical coupling as $\lambda_{c}^{f}\propto\sqrt{K}$, pointing out that the length of the hysteresis cycle increases with $K$, since $\lambda_{c}^{b}$ remains bounded ($\lambda_{c}^{b}\rightarrow 1$). This effect is numerically illustrated in the bottom panel of Fig. 22 where the forward continuation branch $r(\lambda)$ is shown for several star graphs with different values of $K$.

Another interesting approach to the analysis of the star network was presented in Vlasov et al. (2015). In particular, Vlasov et al. noticed that the Kuramoto equations for the star, transformed as in Eq. (42), allow a direct application of the approach introduced in Watanabe and Strogatz (1993, 1994). This way, the set of $K$ coupled equations can be studied by analyzing a single evolution equation for a complex number $z(t)$ whose modulus $|z|$ is the order parameter $r$

[TABLE]

By looking for the steady state solutions of Eq. (47) the authors obtain three such steady state solutions, namely

[TABLE]

Notice that the moduli of these three solutions, $|z^{b}|$, $|z^{f}|$ and $|z^{u}|$ correspond to the three possible values for the order parameter $r$. After performing the stability analysis of the three solutions the authors show that only $z^{b}$ and $z^{f}$ are stable, so that $|z^{b}|$ corresponds to the backward (upper) synchronization branch and $|z^{f}|$ to the forward (lower) one. Moreover, the existence bounds of these solutions ($\lambda^{b}_{c}$ and $\lambda^{f}$) are exactly the expressions in Eq. (46) and Eq. (44) respectively, thus recovering the limits of the hysteresis cycle previously reported.

III.4.3 Effect of structural properties

Up to now we have reported the results of the CDF Kuramoto model in ER and SF networks, i.e. focusing on the effect of degree heterogeneity on the appearance of ES. However, real networks display features like degree-degree correlations between adjacent nodes and modular architectures that may have strong influence in the dynamical processes at work. For this reason, additional research has addressed the role that these two structural properties have on the onset of ES.

Assortativity.- The role of degree assortativity has been studied in several works Liu et al. (2013); Li et al. 2013b ; Sendiña-Nadal et al. (2015). In particular, in Sendiña-Nadal et al. (2015), the authors made an exhaustive study by constructing scale-free networks with tunable degree-degree correlations. The degree-degree correlations were characterized using a Pearson correlation coefficient between the degrees of adjacent nodes as introduced in Newman 2003a

[TABLE]

The networks were produced by a targeted randomization of graphs generated both from preferential attachment and configurational models. The randomization procedure rewires links of the original networks in order to achieve the desired assortativity value $\alpha$ and, importantly, preserves the degree distribution of the original graphs. This way, the authors generate ensembles of networks each one characterized by two parameters: $\gamma$ (the exponent of the power-law degree distribution) and the assortativity $\alpha$.

The numerical simulation of the CDF Kuramoto model unfolding on these networks produces different synchronization diagrams $r(\lambda)$ depending on the values $\gamma$ and $\alpha$. In Fig. 23 (a)-(b) the total area of the hysteresis cycle is reported as a function of $\alpha$ for different scale-free networks. The most important effect is that the assortative networks ($\alpha>0$) display, for moderate values of $\alpha$, much larger hysteresis cycles that their dissasortative counterparts. Interestingly, the increase of degree heterogeneity also magnifies the irreversible nature of ES. Apart from the effect on the hysteresis cycle, Fig. 23 (c)-(d) reveal that the critical couplings, $\lambda_{c}^{f}$ and $\lambda_{c}^{b}$, are also affected by the presence of degree correlations.

Modular networks.- Reported in Li et al. 2013a is the systematic study of ES transition in networks with communities. In this work, the authors study a model network composed of several densely connected modules in which the mixing between modules (the fraction of links connecting nodes belonging to different communities) can be tuned continuously. Their numerical simulations point out that when the modularity of the network is large (so that the mixing between modules is small) the onset of synchronization is enhanced or, equivalently, the value of $\lambda_{c}^{f}$ decreases. Interestingly, to explore in more detail the effects of mixing between different modules, the authors of Li et al. 2013a study a system composed of two star graphs in which the two hubs share some fraction of the leaves surrounding them. The authors observe a double sharp synchronization transition whose position and synchronization gaps can be tuned by controlling the mixing (i.e., the fraction of common leaves) between the two stars.

III.4.4 More general CDF scenarios

The investigations reported on above have assumed a linear correlation between degree and frequency, $\omega_{i}=k_{i}$. However, different correlation functions $\omega_{i}=f(k_{i})$ can be used. In Liu et al. (2013) the authors explore a function of the form: $\omega_{i}=\alpha k_{i}^{\beta}$, so that correlations different from linear can be analyzed. The authors focused on the exponents in the range $\beta\in[-1,1]$, thus covering the case of negative correlations between frequencies and degrees. At first sight, the results show that for $\beta<0$ the abrupt transition is suppressed, becoming second order, and the onset of synchronization is enhanced (see Fig. 24(a)-(b)). Interestingly, the microscopic view of the synchronization transition for $\beta<0$ unveils a very different scenario than that of the the typical second-order transitions found in scale-free networks. Instead of finding a sequential entrainment of oscillators into the synchronized set moving at the average frequency (similar to that shown in Fig. 20(b)), the authors observe for $\beta<0$ what they term an hierarchical synchronization transition.

This hierarchical synchronization transition is illustrated in Fig. 24(c)-(d) for $\beta=-0.1$ and $-0.8$ by monitoring the evolution with the coupling $\lambda$ of the effective frequencies of nodes, Eq. (32), and the average effective frequencies of nodes of the same degree, Eq. (33). These panels show that the oscillators do not evolve monotonously from their natural frequency to the average frequency $\Omega$. On the contrary, the synchronization cluster is governed by those nodes with the largest degrees (and lower natural frequencies). Thus it is the nodes of lower degree that, as $\lambda$ increases, adapt their effective frequencies to join this cluster. As $\lambda$ increases, the average frequency of this synchronized cluster increases, finally reaching the frequency $\Omega$.

A further generalization of the CDF function is analyzed in Skardal et al. (2013). In this case the authors extend the frequency-degree correlation proposed in Gómez-Gardeñes et al. 2011a to $\omega_{i}=\pm\alpha k_{i}^{\beta}$, i.e. they consider the modulus of the natural frequency of a node $i$ as in Liu et al. (2013), but the sign of this frequency is randomly assigned. This way, the probability of finding a node of degree $k$ and frequency $\omega$ becomes

[TABLE]

The purpose of this CDF scheme is to have a symmetric frequency distribution, $g(\omega)=g(-\omega)$, as in the original Kuramoto model, which is a very convenient situation for most theoretical studies.

Under the correlation scheme proposed in Skardal et al. (2013) no explosive transition to synchronization is observed. However, the authors report a two-stage synchronization transition. First, after crossing a first critical coupling, $\lambda_{1}$, a regime of partial synchronization, containing two fully-synhcronized clusters, appears. This is explained as a consequence of the bi-modal frequency distribution $g(\omega)$ obtained from the uni-modal degree distribution. This way, in this regime (termed Standing Wave) the two clusters pace at opposite frequencies, e.g. for $\beta=1$: $\Omega^{+}\simeq\alpha\langle k\rangle$ and $\Omega^{-}\simeq-\alpha\langle k\rangle$. When the coupling increases further, for $\lambda>\lambda_{2}$, the whole network synchronizes. Interestingly, for $\beta=1$ the synchronized regime that appears for $\lambda>\lambda_{2}=2\alpha$ shows that all the oscillators have become phase-locked.

III.4.5 Noisy degree-frequency correlations

The former results are obtained considering a perfect match between degrees and natural frequencies. However, in most real systems this is not the situation and it is interesting to study the CDF model when $\omega_{i}=k_{i}+\xi_{i}$, where $\xi$ is a random small perturbation, $\xi\in(-\epsilon,\epsilon)$. Several works have addressed this issue Zou et al. (2014); Skardal and Arenas (2014).

The most remarkable result in this direction was the stimulation of explosive transitions in networks that, due to their degree homogeneity, do not show this phenomenon in the CDF model. This result, presented in Skardal and Arenas (2014), was shown for one real network, the C. elegans neural network Achacoso and Yamamoto (1991), and two synthetic networks: SF graphs with $\gamma>3$, and networks with a degree distribution given by a stretched exponential (SE), $P(k)=k^{\beta-1}{\mbox{exp}}[-(k/\mu)^{\beta}]$. These graphs show a mild heterogeneity and, as it is proven for SF networks with $\gamma>3$, one should not expect ES to happen. In Fig. 25 the phase diagram for SF ($\gamma=3.5$) and SE networks is shown. In both cases, when the amount of noise in the degree-frequency matching is larger than some critical value $\epsilon_{c}$ ES appears.

III.4.6 CDF with time delay and frustration

Another avenue for testing the robustness of the CDF framework is to consider the effects of time-delay in the coupling between neighbors. Time-delay has been already considered in the original Yeung and Strogatz (1999) and networked Eguíluz et al. (2011) Kuramoto model. Here the question to address is to what extent ES can persist under the presence of time delay.

In Peron and Rodrigues 2012b the authors analyze the following equations

[TABLE]

that for $\tau=0$ recover the original CDF Kuramoto model, Eq. (31). Interestingly, the authors observed a non-trivial dependence of the synchronization transition with the time delay. In particular, they implemented Eq. (50) in a SF network for $\tau=0$, 0.1, 0.3, 0.5, 1.0, 1.5 and 2.0 and discover that only the cases $\tau=0$ and 1.0 lead to an explosive transition. In addition, while the rest of the $\tau$ values show a smooth transition the observed scaling of the order parameter $r$ is rather different.

To gain insight about this anomalous behavior the authors analyze the star-graph. By performing a similar analysis as that in Gómez-Gardeñes et al. 2011a for the existence of phase-locked solutions, the authors find that the expression of the critical coupling $\lambda_{c}^{b}$ needed for the existence of the synchronized state

[TABLE]

where $\omega$ is the frequency of leaf-nodes and $\Omega(\tau)$ the average synchronization frequency. For $\tau=0$ this expression recovers that in Eq. (44). However, when $\tau>0$ the value depends strongly on $\Omega(\tau)$. By using $\omega=1$ for the leaves and $\omega_{h}=K$ for the hub, as in Gómez-Gardeñes et al. 2011a , the authors derive a self-consistent equation for $\Omega$

[TABLE]

From this equation it is clear that the non-monotonous behavior of $\Omega(\tau)$ (and thus the non-trivial behavior of $\lambda^{b}_{c}(\tau)$) is due to the periodicity of the sine function.

Related to these results, the Kuramoto model with frustration (referred to as the Sakaguchi-Kuramoto (SK) model Sakaguchi and Kuramoto (1986)), has been analyzed with the CDF framework Huang et al. (2016); Chen et al. (2017); Kundu et al. (2017). The frustration term modifies Eq. (31) to become

[TABLE]

where $\alpha$ is the frustration parameter. It is important to recall that the SK model was previously found to capture the dynamics of weakly coupled Josephson junctions arrays Wiesenfeld et al. (1996, 1998).

Under synchronization conditions, it is easy to show that Eq. (53) is equivalent to Eq. (50): when the oscillators pace at a common synchronization frequency $\Omega$ we can write $\theta_{j}(t-\tau)=\theta_{j}(t)-\Omega\tau$ and therefore, by comparing Eqs. (50) and (53), we obtain $\alpha=\Omega\tau$. Equation (53) has been analyzed in the star graph in Huang et al. (2016); Chen et al. (2017). In particular, by inspecting the synchronization conditions, in Huang et al. (2016) the authors obtain expressions for both $\lambda_{c}^{b}$ and $\Omega(\alpha)$ similar to those in Eqs. (51)-(52) by changing $\Omega\tau$ by $\alpha$. The study of the SK model in complex networks was tackled in Kundu et al. (2017) finding that, in SF networks, the increase of the frustration $\alpha$ enhances the onset of synchronization at the expense of transforming the ES transition into the usual second order one.

III.5 The correlated frequency-coupling (CFC) framework

The ES transition obtained through the CDF framework motivated the question of whether other mechanisms correlating dynamical and structural aspects of oscillators can lead to same phenomenon. In this regard, Zhang et al. (2013) propose an alternative approach by correlating the natural frequency of an oscillator with the strength of the coupling with its neighbors. This framework, termed here as correlated frequency-coupling (CFC), transforms the Kuramoto model into

[TABLE]

This framework allows us to get rid off the constraint that $g(\omega)=P(k)$ (as required in the CDF) so that one can tackle the study of different frequency distributions $g(\omega)$ using the same underlying network. Why should the CFC framework transform the usual second-order synchronization transition into an explosive one? If we focus on unimodal frequency distributions with mean $\Omega=0$, the effective coupling $\lambda|\omega_{i}|/k_{i}$ for those nodes close to the mean frequency $\Omega$ is very low, compared to that of the nodes which have $\omega_{i}$ that lies on the tails of the distribution $g(\omega)$. Therefore, the CFC mechanism punishes those oscillators that, in the usual Kuramoto model, are the first to get synchronized and, on the other hand, favors those that become phase-locked the last in the traditional setup. In a nutshell, the CFC homogenizes the ability of nodes to become synchronized, thus favoring the sudden emergence of synchronized components.

III.5.1 Independence of network properties

The most remarkable result obtained under this framework Zhang et al. (2013) is that, for a variety of symmetric and unimodal $g(\omega)$ (centered in $\Omega=0$), Eq. (54) displays ES transitions independent of the structural patterns of the underlying network, even for the all-to-all coupling benchmark. In Fig. 26 the forward and backward synchronization diagrams $r(\lambda)$ corresponding to three different network topologies are shown, namely: fully connected network (a)-(b), ER (c) and SF (d) graphs. For the fully connected case, two distributions $g(\omega)$ are considered, Lorentzian and Gaussian, while for the ER and SF graphs the plots refer to Lorentzian distributions. These four synchronization diagrams show the explosive and irreversible nature of the synchronization transition, thus confirming the homogenizing effects of the CFC mechanism.

To shed some theoretical light on the transition the authors in Zhang et al. (2013) propose to follow the approach introduced in Ichinomiya (2004). In a similar fashion as Eq. (34) in the CDF framework, once starts by casting the set of evolution equations for the nodes, Eq. (54) as

[TABLE]

where, as in Eq. (34), $\rho(k;\theta,t)$ is the density of nodes with degree $k$ and phase $\theta$ at time $t$. Here we have neglected degree-degree correlations and the time fluctuations due to finite size effects. By writing the Kuramoto order parameter in terms of $\rho(k;\theta,t)$ Eq. (55) can be rewritten as

[TABLE]

Note that this equation is completely independent of the degree distribution $P(k)$ of the underlying network. This generality comes from the normalization by $k_{i}$ of the coupling between units [see Eq. (54)] as discussed previously in Sec. III.2. At this point it is convenient to work in a reference frame rotating with the average frequency: $\phi(t)=\theta(t)-\Phi(t)=\theta(t)-\Phi_{0}-\Omega t$, so that two phase-locked, $\dot{\phi}=0$, solutions are possible

[TABLE]

These two solutions reveal a novel feature of the CFC framework: the synchronized regime is described as two coexisting clusters of fully phase-locked units. From the form of the two solutions it becomes clear that these two clusters progressively approach each other as $\lambda$ increases. This microscopic picture was confirmed numerically in Zhang et al. (2013) and is shown in Fig. 27.

From the set of solutions in Eqs. (57)-(58) it is possible to evaluate the critical properties of the synchronization transition by expressing the Kuramoto order parameter $r$ in the synchronized regime as the sum of the contributions of the two sets of phase-locked oscillators

[TABLE]

This expression reveals that the phase-locked regime exists when $\lambda>2$, so that the limit of the backward (upper) branch of the synchronization diagram is $\lambda_{c}^{b}=2$. Additionally, the value of $r$ at this point, $r(\lambda_{c}^{b})=\sqrt{0.5}$, indicates the explosive nature of the synchronization transition. This theoretical prediction is, in principle, valid for any network topology. However, for complex networks with small $\langle k\rangle$ the assumptions made to derive Eq. (56) are too strong and, as observed in Fig. 26(c)-(d), fail to predict the correct value for $\lambda_{c}^{b}$.

The analysis of the CFC framework was tackled in several subsequent studies Zhang et al. (2014); Hu et al. (2014); Zhou et al. (2015); Bi et al. (2017); Xiao et al. (2017); Lofti et al. (2018), see Sec. IV for the analysis of the results in Zhang et al. (2014). In Hu et al. (2014) the authors provide with the exact solution of the CFC model in the fully connected graph. Following a similar analysis as in Zhang et al. (2013) the authors reproduce the backward transition finding the two phase-locked solutions $\phi_{1}^{*}$ and $\phi_{2}^{*}$ as in Eqs. (57)-(58) and the critical coupling $\lambda_{c}^{b}=2$. In addition, they are able to solve the forward transition, unveiling how the incoherent state $r\simeq 0$ turns unstable at the value $\lambda_{c}^{f}$. In particular, for even frequency distributions $g(\omega)$, the expression to calculate $\lambda_{c}^{f}$ is the following transcendental equation

[TABLE]

with $\lambda_{c}^{f}$ being the value of $\lambda$ that satisfies Eq. (60) when the real part of the complex number $z$, ${\mbox{Re}}(z)$, changes from negative to positive, i.e. when the incoherent state loses its stability. Importantly, in contrast with $\lambda_{c}^{f}$, the value of the critical coupling for the forward branch depends on the frequency distribution $g(\omega)$ at work.

III.5.2 The Sakaguchi-Kuramoto model under the CFC framework

The SK model was also tackled under the CFC framework in Xiao et al. (2017). In this work the effects of the frustration $\alpha$ over the synchronization transition is shown to strongly depend on the nature of $g(\omega)$. Remarkably, for both triangular and uniform distributions of $g(\omega)$, the increase of frustration implies the inhibition of ES, pointing out the same result observed for the CDF framework in Kundu et al. (2017). However, in the CFC scheme with frustration, the use of a Lorentzian frequency distribution delays considerably the onset of ES and, more importantly, increases the bistability region.

Related to the SK model, the Kuramoto model in the CFC framework with time delays has been recently tackled in Wu et al. (2018). At variance, with the CDF approach, in this case ES is not altered by the addition of time delays (even in the case when they are distributed across the units). Moreover the authors in Wu et al. (2018) show analytically and numerically that time delay can shift the critical coupling $\lambda_{c}{f}$, thus enlarging the bistability region.

III.5.3 Loss of rotational invariance

The presence of the absolute value, $|\omega_{i}|$ in the effective coupling between oscillators in Eq. (54) has very important consequences for the onset of ES. In the seminal article Zhang et al. (2013) the authors noted that replacing $|\omega_{i}|$ by $\omega_{i}$ in Eq. (54) may lead to second order transitions in some situations, thus losing the generality of the results described above. In particular, the transition in Fig. 26(b) changes to second order when the aforementioned replacement in the effective coupling is done. As a conclusion, the presence of the absolute value $|\omega_{i}|$ is key for the CFC framework.

As already mentioned in Sec. III.1, one of the key properties of the usual Kuramoto model, the invariance the equations under rotational transformations, $\theta^{{}^{\prime}}_{i}=\theta_{i}-\omega t$ and $\omega_{i}^{{}^{\prime}}=\omega_{i}-\omega$, is lost due to the inclusion of the absolute value $|\omega_{i}|$ in the effective coupling. Due to this, one should expect non-trivial effects when changing the center of the (even and unimodal) frequency distribution $g(\omega)$ from $\Omega=0$ to $\Omega\neq 0$. This was the object of study in Zhou et al. (2015). The main finding of the authors is that, under the CFC scheme, explosive transitions may turn second order when the center of $g(\omega)$ is translated to non-zero values of $\Omega$.

Interestingly, although the transition becomes second-order when $\Omega$ is large enough, the authors find that the microscopic synchronization patterns are very different from those expected in the traditional Kuramoto model. For small $\Omega$ (when the transition is still explosive) the authors find the two phase-locking clusters reported in Zhang et al. (2013), Eqs. (57)-(58). However, for the second order transition they find (for the forward and backward continuation) two radically different synchronization patterns. Namely, for a small (but nonzero) degree of synchronization those oscillators with small natural frequencies become synchronized forming a cluster of phase-locked oscillators, (i.e., a picture similar to that of the original Kuramoto model). Surprisingly, as $\lambda$ and $r$ increase, two phase-locking clusters show up involving those nodes with large natural frequencies, a pattern similar to that found for the explosive transition. The existence of these two regimes can be derived by writing Eq. (56) in a new reference frame rotating with $\Omega$, so that $\theta^{{}^{\prime}}=\theta-\Omega t$ and $\omega^{{}^{\prime}}=\omega-\Omega$. This way one obtains

[TABLE]

implying that the phase-locked solutions must satisfy

[TABLE]

Note that the function in the l.h.s has the following properties: (i) it diverges at $\omega^{\prime}=-\Omega$, it has a single minimum at $\omega^{\prime}=0$ and it tends asymptotically to $1$ when $\omega^{\prime}\rightarrow\pm\infty$. The form of this function implies that: (i) for small degree of synchronization $\lambda r<1$ only oscillators with $\omega^{\prime}\simeq 0$ ($\omega\simeq\Omega$) will satisfy the condition (62), while (ii) when $\lambda r>1$ the condition is satisfied by those oscillators with positive and negatives frequencies far away from $\omega^{\prime}=-\Omega$. The first condition captures the existence of a single synchronized cluster, while the second one shows the conditions for the emergence of two different clustered of phase locked units.

III.5.4 Including frequency mismatch

An extension of the CFC method was introduced in Leyva et al. 2013b by weighting the coupling between adjacent oscillators by the absolute value of their frequency mismatch. In particular, the authors considered the following equation

[TABLE]

Under this framework, and for $\alpha>0$, those connected units with close frequencies, that are prone to synchronize in the traditional Kuramoto model ($\alpha=0$ within this framework), become weakly coupled. In contrast, the coupling between neighbors with very different frequencies is strengthened by the new term.

The authors studied the emergence of explosive phenomena for different values of $\alpha$ both for homogeneous and heterogeneous networks. For homogeneous, ER and random regular graphs, their results point out that ES naturally shows up for large enough values of $\alpha$ and for a variety of frequency distributions $g(\omega)$. However, by deforming ER into heterogeneous SF networks following the method introduced in Gómez-Gardeñes and Moreno (2006), they observed that the explosive transition is lost progressively as the SF limit is approached and heterogeneity is increased.

III.5.5 Effective Synchronization centrality

To shed more light on how the CFC framework induces ES, in Navas et al. (2015) a method for analyzing how the modification of the coupling between units affects the onset of synchronization clusters is introduced. This method relies on the computation of the eigenvector centrality of each node by considering the modified adjacency matrix corresponding to the effective coupling at work. For instance, in the case of Eq. (63) this effective matrix is $\tilde{A}_{ij}=A_{ij}|\omega_{i}-\omega_{j}|^{\alpha}$, while for Eq. (54) it reads $\tilde{A}_{ij}=A_{ij}|\omega_{i}|/k_{i}$. Once the effective matrix ${\bf{\tilde{A}}}$ is constructed, a subsequent transformation is performed to derive a new matrix as

[TABLE]

with $\Delta\omega_{max}$ being the maximum value of frequency mismatch between connected nodes in the system. The components of the eigenvector of ${\bf C}$ corresponding to the largest eigenvalue, $\vec{v}_{max}$, reveal those nodes that act as synchronization seeds, i.e. those forming the first microscopic synchronization clusters. Thus, the ranking obtained from $\vec{v}_{max}$ of matrix ${\bf C}$ can be seen as a measure of synchronization centrality for nodes.

With this method the authors compare the centrality of nodes in the CFC framework with that of the network with the original adjaceny matrix, i.e. computing $\vec{v}_{max}$ from a matrix ${\bf C}$ derived from ${\bf A}$. The main finding is that, for those networks in which the CFC scheme leads to an ES transition, the CFC acts by homogenizing the synchronization centrality across nodes. This way, the CFC framework effectively suppress those loci originally provided by the network for the creation of synchronization seeds.

III.6 The frequency gap conditioned (FGC) framework

The rationale behind the last approach to the CFC framework, Eq. (63), is to implement an effective decoupling of those connected nodes with similar natural frequencies. This idea is at the core of the next method to be discussed for obtaining ES transitions in networks: the frequency gap conditioned (FGC) approach. However, instead of acting through the evolution equations, as in the cases of the CDF and CFC frameworks, in Leyva et al. 2013a the authors propose to act directly on the design of the network. Namely, given a network topology encoded in an adjacency matrix ${\bf A}$ and a given frequency distribution $g(\omega)$, the aim is to assign natural frequencies to each node in such a way that there is a prescribed mismatch between them given by

[TABLE]

Note that we have set the term $A_{ij}$ so that the mismatch only applies to neighboring nodes.

Obviously, satisfying condition (65) for any given network topology and frequency distribution $g(\omega)$ is hard to achieve and one would face the problem of frustration. To avoid this problem, the authors propose a method to construct networks satisfying Eq. (65). The procedure is as follows: (i) Start from an isolated set of $N$ nodes, (ii) assign randomly the frequencies of each node following the prescribed distribution $g(\omega)$, (iii) choose at random one of the $N(N-1)/2$ pairs of nodes and create a link provided they fulfill condition (65), and (iv) repeat step (iii) until a number of $L$ links have been created. The networks constructed with this procedure display ER-like interaction backbones with $\langle k\rangle=2L/N$.

The authors analyze in detail the role of the frequency mismatch $\gamma$ for a wide variety of frequency distributions $g(\omega)$. In particular, Fig. 28 shows the results corresponding to an homogeneous distribution in the interval $\omega\in[0,1]$ for several values of $\gamma$ (panel a) and $\langle k\rangle$ (panel b). From these plots the authors observe that when the frequency mismatch exceeds some critical value, $\gamma>\gamma_{c}$, ES is observed. The value of $\gamma_{c}$ is reported together with the hysteresis region in panels (e) and (f) for networks with $\langle k\rangle=20$ and $60$ respectively.

Interestingly, when $\gamma>\gamma_{c}$ the resulting networks automatically show degree-frequency correlations, similar to those imposed ad hoc in the CFD scheme. In the top panel of Fig. 28(d) the V-shaped curve reveals these frequency-degree correlations that are necessary to cause the explosive transition and that for $\gamma=0$ (see top panel of Fig. 28(c)) do not show up. In Zhu et al. (2013) the authors explore a very similar FGC framework, making an exhaustive study of how these V-shaped frequency-degree correlations appear exactly at $\gamma_{c}$ and, as a byproduct, they are accompanied by the onset of negative degree-degree correlations in the interaction network.

III.7 The adaptive framework

Up to now, the different frameworks (CDF, CFC and FGC) favoring the onset of ES have focused on coupling the (static) natural frequencies of the oscillators with the structural attributes of the interaction network, either at the level of nodes (CDF) or links (CFC and FGC). A different avenue was proposed in Zhang et al. (2015) by modifying the Kuramoto model to incorporate a new dynamical variable modulating the coupling of a unit $i$ with its neighbors. In particular, the proposed model reads

[TABLE]

where $\alpha_{i}(t)$ can take two different values. Namely, a fraction $f$ of randomly chosen nodes have

[TABLE]

i.e. the coupling is weighted by the instantaneous local field $r_{i}(t)$, defined as in Eq. (24) and normalized by the degree $k_{i}$ so that $r_{i}(t)\in[0,1]$. The remaining $(1-f)N$ units have $\alpha_{i}=1$ as in the original Kuramoto model. Note that for $f=0$ the usual Kuramoto model is recovered.

The formulation of this framework seemingly contradicts the rationale behind the former schemes leading to ES scenarios. Namely, the former three frameworks established different suppressive rules for the formation of small clusters of phase-locked units that serve as the microscopic seeds leading to a large synchronized component. On the contrary, Eq. (66) apparently dictates that those $fN$ oscillators with $\alpha_{i}(t)=r_{i}(t)$ tend to reinforce their connections when the degree of local synchrony is large, thus boosting the creation of the aforementioned microscopic synchronization clusters. However, this reinforcing effect has as its counterpart the null coupling, regardless the value of $\lambda$, when a node is surrounded by dynamically incoherent units. Therefore, the reinforcement of synchronized neighborhoods is compensated by their isolation from the rest of the network.

The numerical findings in Zhang et al. (2015) confirm the above reasoning. In Fig. 29 the forward and backward synchronization diagrams are reported for an ER network when $f=1$. Interestingly, the insets reveal the emergence of ES as a function of $f$ by reporting the average length of the hysteresis cycle $\langle d\rangle$ and the evolution of the two critical couplings $\lambda_{c}^{f}$ and $\lambda_{c}^{b}$ as $f$ varies from [math] to $1$. From these plots it is clear that there exists a critical value $f_{c}$ so that for $f>f_{c}$ the second order transition becomes instead abrupt.

The theoretical explanation for the ES transition in the particular case $f=1$ was presented in Zhang et al. (2015). However, in Danziger et al. (2016) the authors develop a theoretical treatment for the full adaptive model. To this aim, the authors write the effective equations for the phases of adaptive oscillators ($\alpha_{i}(t)=r_{i}(t)$), termed as class I, and of Kuramoto ones ($\alpha_{i}=1$), termed as class II

[TABLE]

where $\phi_{i}$ is the average phase of the neighbors of node $i$. The authors then make a mean-field consideration and approximate the local degree of synchronization to the global one, $r_{i}\simeq r$, and the average phase of the neighbors of $i$ to the global average phase, $\phi_{i}\simeq\Phi$. Note that this approximation is valid when dealing with homogeneous networks having relatively large $\langle k\rangle$. Under these assumptions and changing to a reference frame rotating with the average frequency $\Omega$ (that can be placed at $\Omega=0$), $\psi_{i}=\theta-\Phi$, the phase-locked solutions of Eqs. (68)-(69) satisfy

[TABLE]

Now, considering that in the new rotating frame $r=\sum_{i}\cos\psi_{i}$ and neglecting the contribution of the drifting oscillators (those not satisfying the former equations) to $r$, it is possible to write the following self-consistent equation for the global degree of synchronization that, in the continuum limit, becomes

[TABLE]

The two integrals must be solved by specifying the frequency and degree distributions at work $g(\omega)$ and $P(k)$.

To recover the example of Fig. 29 shown in Zhang et al. (2015) the authors of Danziger et al. (2016) solved the associated integrals numerically finding a surprising result: for any nonzero value of $f$ synchronization becomes explosive. This theoretical result contrasts with the results shown in the insets of Fig. 29. However, the authors note that for small values of $f$ the explosive regime is very small so that, for finite networks is practically undetectable.

The adaptive framework shows its efficiency in causing ES transitions in homogenous (ER) networks. However, its performance in SF networks is still under study. In fact, the architecture of SF networks poses the following question: given the topology constraints imposed by degree heterogeneity and the associated presence of hubs, what is the possibility of achieving the dynamical isolation of small sets of synchronized clusters?

The adaptive framework has been recently reformulated by considering link adaptivity. In this case the connection of the adaptive scheme with the suppresive rule of the static framework is clear. In Avalos-Gaytán et al. (2018) the authors implement an anti-Hebbian rule that enhances the links connecting dynamical units poorly synchronized and weakens the links between coherent oscillators. This framework implements a suppressive rule aimed at delaying the synchronization onset that, as a byproduct, leads to an ES scenario. In addition, the dynamical evolution of link weights naturally evolves towards a CDF architecture, thus explaining the observation of ES transitions.

III.8 Other paradigmatic dynamical models

Most of the studies on ES have been restricted to the Kuramoto model. However, some studies have extended the analysis of ES and the different frameworks addressed above to other well-known dynamical models. Here we report the most relevant ones.

III.8.1 Second-order Kuramoto model

The study of the usual Kuramoto model with inertia has been addressed in several studies Tanaka et al. 1997a ; Tanaka et al. 1997b ; Acebron and Spigler (1998); Olmi et al. (2014). Its implementation in complex networks is also of interest due to its important application to the emergence of synchronous states in power grids Rohden et al. (2012). The usual form of this model reads as follows

[TABLE]

As briefly mentioned in Sec. I.1 the synchronization transition of the Kuramoto model with inertia is one of the few examples of first-order transitions obtained without the addition of terms coupling structural and dynamical features. Therefore, the relevant question to address is how the different frameworks leading to explosive phenomena in the networked Kuramoto model affect the second-order transition.

This question was tackled in the context of the CDF framework in Ji et al. (2013, 2014); Peron et al. (2015); Cao et al. (2018) by defining the natural frequency in Eq. (73) as $\omega_{i}=D(k_{i}-\langle k\rangle)$. The first numerical observations pointed out that, at variance with the CDF Kuramoto model, ES is changed by a cascade of sudden synchronizations of different group of nodes. The authors termed this kind of transition as cluster synchronization and the main novelty is that the sets of nodes synchronizing together are those having the same number of neighbors, i.e., the degree classes. This cluster synchronization is revealed by the sudden changes in $r$ shown in the forward branch of the synchronization diagram.

The cluster synchronization phenomenon is analytically explained in Ji et al. (2014). As for the CDF Kuramoto model, see Eq. (34), the authors consider the continuum limit approach by Ichinomiya (2004). Using the density of oscillators with degree $k$ that have phase $\theta$ at time $t$, $\rho(\theta,t;k)$, and assuming degree uncorrelated networks, Eq. (73) can be expressed as

[TABLE]

where $r$ is the Kuramoto order parameter and, as usual, $\phi$ is the average phase of the ensemble of oscillators.

After changing to a rotating frame that oscillates with the same pace as $\phi$, with $\psi=\theta-\phi$, the non-dimensional version of Eq. (74), reads

[TABLE]

where $\tau=\sqrt{\lambda kr}t$, $\beta=\alpha/\sqrt{\lambda kr}$, $I=D[k-\langle k\rangle-C(\lambda r)](\lambda kr)$ and $C(\lambda r)=(\ddot{\phi}+\alpha\dot{\phi})/D$. This equation is easily recognized as that of a damped pendulum driven by a constant torque $I$. This system is known to display hysteresis when, for instance, at a fixed damping value $\beta$, one increases and then decreases the value of $I$ Strogatz (1994). In particular, the transition from the stable fixed point (here the synchronized state) to the limit cycle (here the incoherent state) and vice-versa do not happen at the same value of $I$.

Following the analogy with the driven pendulum, it is possible to understand the roots behind the abrupt synchronization transition of those group of nodes of the same degree, i.e., the cluster synchronization transition. Recall that each degree class $k$ obeys a different Eq. (75), since $\tau=\tau(k)$, $\beta=\beta(k)$ and $I=I(k)$. To illustrate the different behavior of each degree class, Fig. 30 shows the phase diagram of the pendulum in the $(\beta,I)$-plane for (a) increasing (forward) and (b) decreasing (backward) $\lambda$ (note that bi-stability appears in the region where the limit cycle occurs for increasing $\lambda$ while the fixed point is the local attractor when decreasing $\lambda$). The trajectories corresponding to two degree classes ($k$=8 and $k=20$) are superposed. From Fig. 30(a) is clear that those nodes with degree $k=8$ reach the synchronized state for a smaller value of $\lambda$ in the forward diagram, while Fig. 30(b) shows that nodes with $k=20$ leave the synchronized state before than those of $k=8$.

The scenario of cluster synchronization, however, is not robust under external perturbations. Interestingly, Cao et al. Cao et al. (2018) have recently shown that when the second-order Kuramoto model with CDF is subjected to stochastic perturbations cluster synchronization may be suppressed. Importantly, in these cases, cluster synchronization is replaced by the usual ES of the first-order Kuramoto model.

III.8.2 Oscillators with amplitude dynamics

We now move to the more realistic case of coupled units with both amplitude and phase as dynamical variables, $u_{i}(t)=\rho_{i}(t){\mbox{e}}^{{\mbox{i}}\theta_{i}(t)}$. A paradigmatic model for studying the connection between synchronization of purely phase oscillators with that of phase and amplitude ones is the Stuart-Landau (SL) model Aranson and Kramer (2002),

[TABLE]

where $\omega_{i}$ represents the natural frequency of the $i$-th oscillator. Interestingly, when parameter $\alpha$ is large (compared to the coupling strength $\lambda$) the amplitudes of the oscillators are steady $\dot{\rho_{i}}=0$ and Eq. (76) becomes that of the Kuramoto model, Eq. (23). Thus, the SL model allows a direct connection with Kuramoto dynamics Gambuzza et al. (2016) by tuning parameter $\alpha$. In this class of models two synchronization order parameters are at work, namely,

[TABLE]

The first, $r^{\theta}$, measures the phase coherence while the second, $r^{u}$, measures the degree of global synchronization including the coherence of amplitudes.

The application of the CDF framework in Eq. (76) was studied in Bergner et al. (2012); Gambuzza et al. (2013) for star graphs and complex networks. However, the match $\omega_{i}=k_{i}$ in Eq. (76) does not produce ES in any of the topologies studied. On the other hand, the authors observed that the CDF scheme produces a novel phenomenon called remote synchronization which refers to phase coherence between units that are neither directly connected nor connected by a finite path of units sharing phase coherence.

At variance with CDF, the implementation of the CFC framework in the SL model gives rise to explosive transitions. However, in this case the abrupt transition refers to the sudden death of oscillations (OD). This result was presented in Bi et al. (2014) where the authors implemented the frequency-coupling correlation in Eq. (76) by changing $\lambda\rightarrow\lambda|\omega_{i}|/N$. The authors focused on fully connected graphs and analyzed different frequency distributions $g(\omega)$: triangular, Lorentzian and homogeneous. In the three cases the numerical simulations showed explosive transitions with bi-stability regions for both $r^{\theta}$ and $r^{u}$. However, the microscopic analysis of these explosive transitions yield some differences in terms of how the effective frequencies of oscillators approach to the mean $\Omega=0$. On the other hand, the common feature of these transitions is that, in the synchronized regime, the oscillators self-organize into two clusters of the same amplitude but with two different phases. In Fig. 31 the three explosive transitions are reported via the diagrams $r^{u}(\lambda)$, the microscopic evolution with increasing $\lambda$ of the effective frequencies $\omega^{\mbox{\scriptsize eff}}_{i}$, Eq. (32), and the difference between the phases of the oscillators and the average phase $\Phi$.

Apart from the SL, other paradigmatic models of coupled oscillators with phase and amplitude dynamics have been analyzed. In Chen et al. (2015), the Complex Ginzburg-Landau Aranson and Kramer (2002) is analyzed in a population of homogeneous ($\omega_{i}=\Omega$) oscillators interacting through a SF network. The authors show that, despite of the dynamical homogeneity of this system, explosive phase coherence shows up in a wide region of parameter space. In Verma et al. (2017) and Zhao et al. (2018) explosive transitions to amplitude death (AD) are found for coupled Van der Pol (VP) oscillators. In particular, in Verma et al. (2017) a set of mean-field coupled VP oscillators is analyzed while Zhao et al. (2018) tackles the networked version. Interestingly, in the latter work it is shown that the backward transition point is independent of the underlying network topology while the average degree has a strong impact on the nature of the explosive (AD or OD) transition found.

III.8.3 Neural models

Following with equations for diffusively coupled oscillators with cubic nonlinearities, we now report the results in Chen et al. 2013a for the FitzHugh-Nagumo (FN) model Murray (1989). In this work the authors consider the system analyzed in Pikovsky and Kurths (1997) of coupled FN oscillators,

[TABLE]

In these equations $\epsilon=0.01$ so that, as usual in relaxation oscillators, we have a fast ($x$) and a slow ($y$) variable. Finally, $\xi_{i}(t)$ is the term of Gaussian noise of unit $i$ satisfying $\langle\xi_{i}(t)\rangle=0$ and $\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2D\delta_{ij}\delta(t^{\prime}-t)$.

The relevant parameter governing the behavior of each unit is $\alpha_{i}$. In particular, for $|\alpha_{i}|>1$ the unit becomes excitable, while $|\alpha_{i}|<1$ leads to an oscillatory behavior. Since the frequency of unit $i$ is a decreasing function of $\alpha_{i}$ the authors in Chen et al. 2013a applied the CDF framework by defining

[TABLE]

where the parameter $\delta$ controls the slope of the dependence of $\alpha_{i}$ with the degree of the node. Given that $\delta>0$, from Eq. (81), is clear that the oscillation frequency of a node $i$ increases with its degree.

By solving numerically Eqs. (79)-(80) the authors study the synchronization diagrams for ER and SF networks for different values of $\delta$. These simulations, see Fig. 32, reveal that only SF networks show ES transitions. However, a careful inspection of the role played by $\delta$ (see Fig. 32(g)) highlights that there is critical value $\delta_{c}\simeq 0.2$ for the onset ES. At variance with most of the models reviewed, for $\delta\gtrsim\delta_{c}$ the length of the bistable region $d$ suddenly takes large values. Moreover, as $\delta$ increases the length of the bi-stable region decreases.

III.8.4 Rossler chaotic oscillators

To complete the discussion, we move to the domain of chaotic oscillators. In Leyva et al. (2012) the authors analyze a system of diffusively coupled Rossler-like units. In particular, the authors analyze the following Rossler system introduced in Pisarchik et al. (2006, 2009),

[TABLE]

[TABLE]

As usual in Rossler units, nonlinearity enters in the $z$ variable of each unit and, for this particular system, it is specified by a piecewise function defined as

[TABLE]

In Leyva et al. (2012) some of the parameters are fixed ($\Gamma=0.05$, $\beta=0.5$, $\epsilon=1$ and $\mu=15$) while $\nu=0.02-\frac{10}{R}$, where $R$ is a tunable quantity that regulates the dynamical state of the system. In particular, for the values specified above, the system operates in the chaotic region when $R\in[55,110]$ Pisarchik et al. (2006). The second tunable parameter is the coupling between units $\lambda$ that, in this case, is applied through the $x$ variable.

Since the natural oscillation frequency of a node $i$ depends linearly on the parameter $\alpha_{i}$, the authors apply the CDF framework by imposing a positive correlation between $\alpha_{i}$ and the degree $k_{i}$ of each unit according to

[TABLE]

where $\alpha=10^{4}$ and, as in Eq. (81), $\delta$ controls the slope of the function $\alpha_{i}(k_{i})$. Note that, at variance with Eq. (81) here $\alpha_{i}$ is not bounded and, therefore, for heterogeneous degree distributions the CDF framework implemented yields a very broad frequency spectrum.

The synchronization diagrams obtained from the numerical solution of Eqs. (82)-(84) are reported in Fig. 33(a) for three SF networks with different exponents $\gamma$. The results for the SF network with $\gamma=3$ point out that, depending on the dynamical region (controlled by $R$), this network shows explosive or continuous transitions. In Fig. 33(b) the synchronization diagram $r(\lambda,R)$ for this network confirms the former result. For $R\gtrsim 88$ (region II) ES of chaotic states of chaos appears, together with a region of bistability.

III.9 Experiments and real-world applications

The variety of mechanisms leading to ES allows for them to be implemented in controlled experimental setups and observed in real-world systems. Here we review the experimental observations of ES reported in the literature.

III.9.1 Implementations in the lab

The numerical observation of ES in the Rossler system allowed the authors of Leyva et al. (2012) to validate experimentally the existence of this phenomenon in a real setup subjected to stochastic fluctuations and parameter mismatches. In particular, they implemented the electronic network shown in Fig. 34.a consisting of $N=6$ piecewise Rössler circuits (labeled as $N1$, $N2$, …, $N6$) having the same internal $R_{\mbox{\scriptsize exp}}$, so to ensure that they operate in an almost identical dynamical (chaotic) regime. These circuits were constructed as in introduced in Pisarchik et al. (2006, 2009) so that they obey the evolution described in Eqs. (82)-(84).

Due to the impossibility of constructing large-scale SF networks as those used in numerical simulations, the circuits were arranged in a star-like configuration, as this topology has been previously shown (Sec. III.4.2) to display ES in the CFD Kuramoto model. To implement the frequency mismatch between the central and peripheral nodes, the authors made a fine tuning of the values of the capacitors so that the central circuit ($N1$) oscillates with a mean frequency of $3,333$ Hz, and those circuits acting as leaves ($N2$, …, $N6$) have frequencies in the range of $2,240$ $\pm$ $200$ Hz. Obviously, the frequencies of the circuits are slightly different due to experimental variability.

The experiment was performed by coupling simultaneously the $6$ chaotic oscillators by means of $5$ digital potentiometers (red nodes) to ensure that the coupling strength $\lambda$ can be tuned simultaneously for all the pairs of connected nodes. This way, for each value of $\lambda$, the output signal from each circuit is stored and the time evolution of the phase of each circuit is reconstructed. These phases are then analyzed to compute the global order parameter $r$ of the system.

In Fig. 34(b) the experimental forward and backward synchronization diagrams, $r(\lambda)$, are reported for $3$ values of the parameter $R_{\mbox{\scriptsize exp}}$. In close analogy with the results for SF networks, Fig. 33(b) shows that the nature of the synchronization transition crucially depends on the value of $R_{\mbox{\scriptsize exp}}$. For $R_{\mbox{\scriptsize exp}}=55$ (just in the limit of the chaotic regime) the transition to synchrony is smooth while the diagram for $R_{\mbox{\scriptsize exp}}=65$ and $70$ shows a sharp increase of the synchronization parameter $r$ together with the hysteresis cycle. This experimental observation was the first one to put the theoretical predictions about ES into practice.

Another interesting experimental set up has been proposed in Kumar et al. (2015) where the authors analyze a star network of mercury beating-heart (MBH) oscillators. A MBH is an electrochemical redox reaction between mercury, iron and chromium discovered by Lippmann (1873). The reaction causes a blob of mercury in water to oscillate. The authors implement the CDF framework in this system and show that, as the coupling strength is increased (forward direction) and decreased (backward direction), the forward and backward curves of the synchronization order parameter do not overlap, thus indicating the existence of a hysteresis (bistability) region, where explosive synchronization emerges. Moreover, it is given evidence that the dynamics of the oscillators can be conveniently and permanently switched from one phase to the other superimposing an instantaneous external signal, if the system operates in the bistable region.

III.9.2 Spontaneous synchronization in real systems

The first reported study on the phenomenon of spontaneous synchronization in networked real systems, was presented in Motter et al. (2013) for the case of the electric power-grid. The authors show that the leading factor for stability is the relation between the specifics of the dynamical units and the network structure. To this aim they represent each dynamical unit using the swing equation Anderson and Fouad (2008), and couple them in a topological network in accordance with real configurations of the power-grid around the world. The object of analysis, however, is not the adjacency matrix of the physical topology of connections but the admittance matrix, defined as the Laplacian matrix for the admittance (a measure of how easily a device will allow a current to flow) between nodes. The stability analysis of the system is performed by assuming small perturbations of the synchronized state, in the same fashion as the master stability function Pecora and Carroll (1998). This analysis results in a condition for the Lyapunov exponents, connecting topology and dynamics. Other previous works have also addressed the synchronization properties of (synthetic) power-grids Lozano et al. (2012); Rohden et al. (2012). The spontaneous synchronization observed in real power-grids is then described as an outcome of the destabilization of the power-grid.

Apart from the previous works, the other interesting field where spontaneous synchrony seems to be important is neuroscience. There is an ongoing consensus in biology Chatterjee et al. (2008) about the role of bistability as a natural framework for biological switches. This bistability is many times reported as an hysteresis cycle controlling the response of the system to external stimuli, with synchronization being a paradigmatic example behind many phenomena. Following these ideas, many abrupt transitions and hypersensitive responses to perturbations observed in neuroscience have been recently revisited from the perspective of ES. Examples include the conscious-unconscious transition Kim et al. (2016, 2017), epileptic seizures in the brain Wang et al. (2017), the sensitive frequency detection of the cochlea Wang et al. (2016), the hypersensitivity of Fibromyalgia patients Lee et al. (2018), and the activity of the choroid plexus controlling mammalian circadian clocks Myung et al. (2018).

In Kim et al. (2016, 2017) the authors hypothesized that the conditions for explosive synchronization in human brain networks would be present in the anesthetized brain just over the threshold of unconsciousness. In particular, network conditions for explosive synchronization have been found in empirically-derived functional brain networks Joiner et al. (2013). Unconsciousness and resting states Deco et al. (2013) seem to be related to a bifurcation point from which the dynamical system confronts multistability that may lead to spontaneous synchrony.

In Wang et al. (2017) the authors investigate the onset of epileptic seizures as the onset of ES. This study reveals that ES may be induced due to a rewiring of the underlying network that increases local clustering.

Also recently, the authors in Wang et al. (2016) proposed a model in which acoustical signal transduction in the cochlea (a spiral-shaped cavity forming the auditory portion of the inner ear) presents an abrupt synchronization transition which provides a plausible explanation for the system’s ability for frequency selectivity.

All of these works are only beginning to reveal the horizon, and more work is still necessary to accurately relate the phenomenon of ES to observations in real systems. One example linking the observation of ES with the presence of one of its related mechanisms is the connection between Fibromyalgia and ES. By analyzing EEG signals from patients with chronic pain, the authors of Lee et al. (2018) revealed a significant positive correlation between the frequency and degree of cortical areas within the alpha frequency, pointing out that ES (under the CDF framework) could be the underlying mechanism behind the hypersensitivity that characterizes this condition. Another example along this same line Myung et al. (2018) has unveiled that the activity in the choroid plexus of mice adjusts the suprachiasmatic nucleus clock, responsible for generating the circadian rhythms. When analyzing the dynamical patterns of gene expression found in the choroid plexus, the authors suggest that the frequency and coupling of individual oscillators in the choroid plexus are positively correlated (CFC framework) and that this mechanism is critical to how synchronization and spatial patterning arise in the choroid plexus.

III.10 ES in a nutshell

To round off this section devoted to ES we summarize and compare the main frameworks leading to this phenomenon. In a similar fashion as for EP (recall Fig. 12) this summary is reported in Fig. 35 where, in the first column, the basic rules of the CFD, CFC, FCG and Adaptive frameworks are shown as they are implemented in the Kuramoto model together with the references of the works introducing each framework.

The toy networks drawn in the second column show how the former rules are implemented. The colors of the nodes represent the natural frequencies of nodes, see legend on the top (adjacent colors represent close values for $\omega$). This way, the assignation of $\omega$ in the CFD model is made according to the degree of each node. In the cases of the CFC, FCG and Adaptive frameworks the frequencies have been randomly assigned. The distribution of frequencies affect, in these cases, the interactions among the nodes. Namely, for the FCG and Adaptive frameworks the backbone network is preserved but the redefinition made for the interaction coupling of a node $i$ with its neighbors turns the original network into a directed one. In addition, the directional links are weakened or reinforced depending on the attributes of the nodes (CFC model) or the dynamical state of the neighborhoods (Adaptive framework). For the FCG case the frequency distribution biases the skeleton of interactions since the existence of a link is restricted to those pairs of nodes with a pre-defined frequency gap (dashed lines).

The third and fourth column of Fig. 35 report the dependence of ES on the frequency and degree distributions. While for the CFD model the degree distribution imposed automatically forces the frequency distribution to be the same, the rest of the ES cases reported have been for a variety of typically even frequency distributions (such as Gaussian, Lorentizian, homogeneous, etc). Note that in the CFC framework the definition of $g(\omega)$ has to be made keeping in mind that, due to the loss of rotational invariance of the equations introduced by the redefinition of the coupling, the distribution has to be centered around $\Omega=0$. The existence of ES is also influenced by the $P(k)$ of the underlying network. For instance, in the CFD model only SF networks with $\gamma<3$ display ES while for the FCG model the $P(k)$ of the network is imposed by the frequency assignment given by $g(\omega)$, giving rise to homogeneous (ER-like) topologies. For the CFC and Adaptive models there is no restriction for $P(k)$ although, for the Adaptive case, the existence ES in SF networks is still to be studied.

Finally, the fifth column in Fig. 35 reports the information horizon needed to implement the different rules leading to ES. For the CFD and CFC this horizon is purely local since only the attributes of single nodes are needed. This horizon expands for the Adaptive framework since the dynamical states of the neighbors of each node $i$ are needed to construct $r_{i}$. For the FCG case the information horizon expands further to the whole network since the construction of the contact network implies the search to all possible pairs of nodes that fulfill the frequency gap condition.

IV Connection between EP and ES

The phenomenology behind EP and ES share similar ingredients, but there are also apparent mechanistic and phenomenological dissimilarities. Both EP and ES are fundamentally based on the implementation of some microscopic rules that avoid the formation of a macroscopic (connected or synchronized) component that leads to the onset of the collective phenomenon at work. In EP these rules are explicit. In ES the rules are implicit and result from the combination of the network topology and the natural frequency distribution of the nodes. Table 3 summarizes some apparent key distinctions between EP and ES. Whereas EP is usually defined as a kinetic process with links between nodes added or removed, ES is based on a dynamical process where nodes are either in sync or not, depending on their phase difference and frequencies. In ES the network dynamics run on a underlying fixed network whereas in EP the underlying structure evolves as nodes across distinct clusters are linked together. Typically, EP is based on a non-local cluster merging process while in ES cluster merging is a local node-to-node based process. EP requires a sufficiently homogeneous initial cluster size distribution whereas in ES the node frequency distribution is usually heterogeneous. Finally, in EP, at $p_{c}$, a gap between small-sized and larger-sized clusters is observed, whereas for ES a necessary condition can be formulated in terms of frequency gaps.

While an exact map between EP and ES is yet to be established, in Zhang et al. (2014) the authors monitor the formation of microscopic synchronized clusters in order to translate the synchronization transition as the coupling strength $\lambda$ varies in terms of a percolation-like process. In particular, the authors apply the technique introduced in Gómez-Gardeñes et al. 2007a to monitor the formation of synchronized links by measuring the local degree of synchronization between two nodes $i$ and $j$ as

[TABLE]

where $T$ denotes the time-interval over which integration is performed. Following this definition, $R_{ij}=1$ when two nodes are fully synchronized, and $R_{ij}=0$ when their dynamics are completely incoherent.

For the Kuramoto model in complex networks, Eq. (23), the computation of $R_{ij}$ as a function of $\lambda$ allows for the study of how synchronized clusters emerge by monitoring the formation of each synchronized link between two nodes, say $i$ and $j$, whenever $R_{ij}>\eta\rightarrow 1$. In Gómez-Gardeñes et al. 2007a the creation of synchronized links in ER and SF networks was monitored using this technique establishing two qualitatively different paths toward synchronization. Namely, as shown in Fig. 36, for ER graphs the transition starts with the formation of separated clusters of synchronized nodes that, at $\lambda=\lambda_{c}$, coalesce forming a macroscopic (in terms of nodes) cluster. For $\lambda>\lambda_{c}$ the clusters grow incorporating links between the already present nodes until full synchronization is reached Gómez-Gardeñes et al. 2011b . On the other hand, for SF networks the synchronization path is governed by a unique cluster initially formed by the hubs. This cluster grows monotonically by incorporating new nodes as $\lambda$ increases.

In Zhang et al. (2014) the authors follow the same procedure but apply it to the CFC Kuramoto model, Eq. (54). This way they implement the forward (backward) continuation in $\lambda$ of the CFC model and measure the correlation matrix ${\bf R}$ for $\lambda$ values around $\lambda_{c}^{f}$ ($\lambda_{c}^{b}$). The results are shown in Fig. 37 for three network topologies with $N=200$ nodes: (top row) fully connected network, (middle row) ER graph, and (bottom row) SF network. In all the cases $g(\omega)$ is Gaussian with $\Omega=0$. The first three panels in each row (starting from the left) correspond to $\lambda$ values before the onset of ES in the forward branch, while the last one shows the matrix ${\bf R}$ immediately after $\lambda_{c}^{f}$. The rows (and columns) are arranged in ascending order of natural frequencies so that nodes with similar $\omega_{i}$ correspond to adjacent rows (columns) in matrix $R$.

The results of this study point out that when $\lambda<\lambda_{c}^{f}$ small synchronization clusters appear. However, these clusters are formed along the diagonal, pointing out that they are formed by nodes of similar frequencies, avoiding the creation of synchronized links connecting nodes of disparate frequencies (far from the diagonal). These clusters remain isolated even at the $\lambda\lesssim\lambda_{c}^{f}$, in contrast to the traditional paths to synchronization observed in complex networks (Fig. 36). Whereas for $\lambda\gtrsim\lambda_{c}^{f}$ all clusters merge together forming a fully synchronized network.

To illustrate how small clusters of synchronized nodes avoid merging together even when $\lambda\lesssim\lambda_{c}^{f}$ the authors in Zhang et al. (2014) derived a suppressive rule for the formation of synchronized links in the CFC Kuramoto model. Consider two oscillators with $\omega$ and $\omega^{\prime}$ as natural frequencies, from Eq. (56) in Sec. III.5.1, the evolution for their phase difference $\Delta\theta=\theta-\theta^{\prime}$ is

[TABLE]

Therefore, these two oscillators will be locked, $\Delta\dot{\theta}=0$, when

[TABLE]

The suppressive rule points out that those oscillators with close frequencies that, in addition, lie on the tails of the distribution $g(\omega)$ (centered at $\Omega=0$) will be those establishing synchronized links while for those around $\Omega=0$ the condition of Eq. (89) is difficult to fulfill. This becomes evident from the absence of synchronized links in the central part of the matrices ${\bf R}$ in Fig. 37 corresponding to $\lambda<\lambda_{c}^{f}$, as these entries correspond to pairs of nodes with close and small absolute value frequencies.

Recently, another study has shed light on the behavior of synchronized clusters at the ES onset Koronovskii et al. (2017). Focusing on the backward transition, they study how the giant synchronized cluster is destroyed at $\lambda_{c}^{b}$. Their results reveal that, both for the CFD and CFC frameworks, the synchronized cluster disappears progressively in time, step by step, keeping the probability distribution of the phases of its constituents invariant during the process. Thus, this finding indicates the self-similarity of the synchronized clusters that are observed at $\lambda_{c}^{b}$.

A mapping between EP and ES phenomena may be enabled by establishing connections between the explicit rules of EP and the implicit rules of ES. For instance, the suppressive rule of Zhang et al. (2014) is one such attempt. It is also worth noting that not all of the explicit rules that lead to EP directly suppress the growth of the largest cluster, for instance the Devil’s staircase rule operates by preferentially merging clusters of similar size Nagler et al. (2012). Likewise, one can initialize a set of clusters to have a “powder keg” and then any edge addition rule, even random addition, leads to EP Friedman and Landsberg (2009); Cho et al. (2010). Similarly, an ES transition may be prepared by a sufficiently broad natural frequency distribution of the nodes, but the sufficiency is dependent on the underlying network topology. All-to-all coupling enhances the conditions for a continuous transition as the fraction of out-of-sync nodes slows the growth of the largest cluster. On the other hand, a network topology with pronounced community structure or degree disassortativity may enable ES. Thus, a framework that unifies explosive transitions appears possible for some but not for all explosive phenomena.

V Other models with explosive behavior

V.1 Explosive congestion in complex networks

Transportation dynamics on networks can be, in general, interpreted as the flow of elements from an origin node to a destination node. When the network is used by simultaneous transportation processes, many elements necessarily travel through the same node or link. This, in combination with the possible physical constraints of nodes and links, such as the capacity to route a packet in finite time, can lead to network congestion, a phenomenon in which the amount of elements in transit on the network grows proportional to time Guimerà et al. 2002b ; Guimerà et al. 2002a ; Zhao et al. (2005). Network congestion is always a function of the network topology and the routing protocol. The most widely used transportation protocol is shortest path routing, a protocol consisting in minimizing the distance (or time) between two places.

The pioneering work by Echenique et al. (2005) explored how the average network performance depends on the ability of the routing protocol to divert traffic across paths other than the shortest ones. To this end, the authors numerically explored a model in which a tunable parameter accounts for the degree of traffic awareness incorporated into a packet’s delivery. Suppose that a node $l$ is holding a packet whose destination is node $j$. Any node $i$ that is a neighbor to $l$, has an effective distance between $i$ and $j$ defined as

[TABLE]

where $d_{i}$ is the minimum number of hops one has to pass by to reach $j$, i.e., the shortest path between $i$ and $j$, and $c_{i}$ is the length of the queue of packets to be delivered by node $i$. Additionally, $h$ is a tunable parameter that accounts for the degree of traffic awareness incorporated in the delivery algorithm. The main finding is that the onset of the congested phase is reminiscent of a first or a second order phase transition depending on whether or not the routing combines a shortest path delivery strategy with traffic aggregation at a local scale. This abrupt transition can be understood as a shifting of the critical point of a second order phase transition, that when attained reveals the large distribution of congestion across many points in the network, see Fig. 38.

The explosive emergence of congestion in scale-free networks was further analyzed by Hu et al. (2007) proposing a model where the node queue length is proportional to the node degree and its delivering ability proportional to the queue length. They found a cycle of hysteresis, indicating that the system is bistable in a certain range of packet density. For different packet densities, the system can self-organize to four different phases: free-flow, saturated, bistable, and jammed. The bistability is related to the existence of a first order phase transition in the flow.

Later on, in De Martino et al. (2009) the authors proposed a minimal model to reveal the mechanisms triggering the emergence of congestion. These mechanisms are reminiscent of a bootstrap percolation process, where a node is active (here, congested) if the number of active neighbors exceeds a given threshold (in this case the length of the queue). A similar finding of the explosive transition to congestion, was exposed in Wang et al. (2009) for finite buffer nodes. The authors propose also control strategies, based on the idea of random packet dropping, to prevent complete congestion. Such control strategies can modify the order of the phase transition from first to second order.

In a more general framework of dynamic information routing Kirst et al. (2016), the authors explore a generic mechanism to route information on top of collective dynamical reference states in complex networks. A dynamical reference state is simply a global solution of the dynamical process (e.g. synchronization). They prove specifically for a Kuramoto-like oscillatory dynamics, how individual unit properties, the network topology and external inputs act together in concert to organize information routing. Interestingly, they found that changing for example the unit properties (in the case of the Kuramoto model, the intrinsic frequency of oscillators), the transition in information routing has a switch-like dependency on the control parameter reminiscent of the explosive transitions above mentioned.

V.2 Explosive epidemics in networks

Epidemic spreading has a long tradition in the physics literature Pastor-Satorras et al. (2015) with a wide interest in the determination of the epidemic threshold: a value of the infectivity that separates two scenarios, non-invasive and invasive regimes, respectively. In particular, the epidemic threshold concept is analogous to the concept of critical point for a phase transition in non-equilibrium systems. This epidemic threshold phase transition has been proven to be continuous and second-order for the vast majority of epidemic dynamical models. However, from the pioneering work by Gross et al. (2006) on adaptive networks, which includes the possibility for individuals to protect themselves by avoiding contacts with infected people, the authors found the presence of a hysteresis loop and first order transitions in the fraction of infected individuals as a function of the rewiring (defined as the probability of redirecting connections to non-infected individuals, $w$), see Fig. 39.

Usually, epidemic models assume that the transition to endemic states in a population can be fully explained in terms of contagions between pairs of individuals. However, several epidemic models are also used to understand social information contagion, and the dynamics of social transmission depend not only on the characteristics of the transmitting and receiving individuals (e.g. on attitude or persuasiveness) but also on the context of the transmission event. In Gómez-Gardeñes et al. (2016), the authors propose a model that accounts for explosive contagion in complex networks. Essentially, the model relies on a definition of the transmission from a transmitter $j$ to a receiver $i$ with synergistic rate given by

[TABLE]

where $\alpha$ is the intrinsic value of the spreading phenomenon in the absence of the context. The function $f(x)$ encodes the synergy of $x$ on the result of the transmissibility, and the argument $n^{h}(i)$ is the number of ignorant/healthy individuals connected with the receiver $i$. For two different choices of $f(x)$, linear and exponential, the authors explore the effect of a multiplicative parameter $\beta$ on $n^{h}(i)$, such that if $\beta>1$ the synergetic effect is positive, while for $\beta<1$ the synergy is destructive. Using a microscopic Markov Chain approach (MMCA) Gómez et al. (2010, 2011), it is possible to analytically describe the evolution of $p_{i}(t)$, which is the probability that a node $i$ becomes infected at time $t$ in a susceptible-infected-susceptible (SIS) epidemic model as follows

[TABLE]

Here $\mu$ is the recovery rate of infected individuals, and $q_{i}(t)$ the probability of being infected by a neighbors, which reads

[TABLE]

where $A_{ij}$ is the network adjacency matrix. The function $\lambda_{j\rightarrow i}$ is thus:

[TABLE]

The phase diagram of the system is depicted in Fig. 40. Both the Monte Carlo simulations and the MMCA theory predict the coexistence of endemic infected and infection-free states and the corresponding hysteresis effect with discontinuous transitions between these regimes. Similar results have been reported in the dynamics of the susceptible-infected-recovered (SIR) model with reinfections Gómez-Gardeñes et al. (2015).

A prescription that also shows the phenomenology of explosive epidemics is found in Liu et al. (2017), with a focus on local synergy, quantified by the interaction of elements that when combined produce a total effect greater than the sum of the individual elements. In the case of epidemic contagion, the authors study a model in which the usual infectivity parameter $\beta$ is rescaled as

[TABLE]

where $\alpha$ is the strength of the synergistic effect, and $m$ the number of infected neighbors. The effect of synergy is an explosive epidemic outbreak characterized by a first order phase transition, see Fig. 41 for details.

A very interesting effect was reported in Böttcher et al. (2015) when the recovery of sick individuals depends on the availability of healing resources that are generated by the healthy population. In this case, the limitation of resources triggers the emergence of an abrupt transition in the epidemic outbreak. Abrupt percolation transitions in complex networks can also be produced by the spreading of cooperative contagions Cai et al. (2015). Let us imagine a set up in which the agents can be infected by a single disease A with probability $p$, and that nodes infected with A can be infected with a different disease B with a raised probability $q$. The net effect of this dynamical process reveals that cooperation can lead to phase transitions of different orders, causing either a mild or an abrupt massive outbreak at the threshold. An abrupt outbreak emerges only if the network structure supports a bottleneck for cascading mutual infections. Other specific situations where epidemics can emerge in an explosive transition are, for example, when the infection spreads and in addition we require that the recovery of infected nodes can only occur when they remain connected to a predefined central node via a path that contains only healthy nodes Böttcher et al. (2016). This particular situation induces the possibility that healthy nodes become encapsulated by infected nodes and thus are absorbed into an infected cluster all at once, thereby promoting abrupt transitions at the epidemic outbreak.

V.3 Explosive phenomena on multilayer complex networks

The explosive phenomena surveyed in the previous sections where grounded in “traditional” complex networks, in that the focus was on individual networks. However, there is an increasing interest in so-called “networks of networks”, accounting for a diversity of distinct network types, node and link types, and interconnection patterns. A useful paradigm for organizing the complexity of many such systems is the framework of multilayer networks Kivelä et al. (2014); Boccaletti et al. (2014), see Fig. 42 for an example. Here we briefly review the findings of abrupt (explosive) transitions that emerge from the structural aspects of a multilayer network’s topology, as well as from the interplay between structure and certain dynamical processes on multilayer networks.

Interdependent and multilayered networks are used to define classes of complex systems where networks are not considered as isolated entities but instead as interacting with one another. More specifically, we can consider a multiplex network, which is a multilayer network in which each distinct network layer implements a different categorical relationship between the entities in that layer, and the same exact set of nodes exists in each network layer as replicas of each other. Note, the notion of interdependent networks is a more general framework where nodes can be different in each network, and the functioning of nodes in one network can depend upon interconnections to nodes in another network, e.g. the electric grid and the cyber-infrastructure used for its control. Many, if not all, real networks are “coupled” with other real networks. Examples can be found in several domains: social networks (e.g., Facebook, Twitter, etc.) are coupled because they share the same actors Szell et al. (2010); multimodal transportation networks are composed of different layers (e.g., bus, subway, etc.) that share the same locations Barthélemy (2011); the functioning of communication and power grid systems might depend one on the other Buldyrev et al. (2010). The majority of phenomena that have been studied on interdependent networks, including percolation Buldyrev et al. (2010); Son et al. (2012), epidemics Saumell-Mendiola et al. (2012), evolutionary games Gómez-Gardeñes et al. (2012) and linear dynamical systems Gómez et al. (2013), have provided results that differ significantly from their counterparts obtained in the case of isolated complex networks De Domenico et al. (2016). Sometimes the difference is striking: for example, while the large-scale connectivity of isolated scale-free networks is robust against random failures of nodes or edges Albert et al. (2000), scale-free interdependent networks can instead be very fragile with respect to the breakdown of large-scale connectivity Buldyrev et al. (2010); Son et al. (2012).

V.3.1 Abrupt structural transitions

The merging of independent networks into an interconnected network-of-networks undergoes a structurally sharp transition as the interconnectivity between the networks increases. Depending on the relative importance of inter- and intra- layer connections, it has been found Radicchi and Arenas (2013); Arenas and Radicchi (2016) that the system can exist in one of two regimes: in one regime, the various layers are structurally decoupled and they act as independent entities; in the other regime, network layers are indistinguishable and the whole system behaves as a single-level network, with an abrupt transition between the two regimes.

Focusing on the analysis of the spectrum of the supra-Laplacian matrix, defined as

[TABLE]

where the blocks present in $\mathcal{L}$ are square symmetric matrices of dimensions $N\times N$. In particular, $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ are the Laplacians of the networks $A$ and $B$, respectively, and $p$ the inter-connectivity weight. The spectrum of the Laplacian of a graph is a fundamental mathematical object for the study of the structural properties of the graph itself. The authors in Radicchi and Arenas (2013) unveiled the origin of the structural changes of the merging of networks in an interconnected system. Note that, for any graph, all eigenvalues of its Laplacian are non-negative numbers. The smallest eigenvalue is always equal to zero and the eigenvector associated to it is trivially a vector whose entries are all identical. The second smallest eigenvalue $\lambda_{2}$, also called the algebraic connectivity, is one of the most significant eigenvalues of the Laplacian. It is strictly larger than zero only if the graph is connected. More importantly, the eigenvector associated to $\lambda_{2}$, which is called the characteristic valuation of a graph, provides even deeper information about its structure. For example, the components of this vector associated to the various nodes of the network are used in spectral clustering algorithms for the bisection of graphs. By studying the behavior of the second smallest eigenvalue of the supra-Laplacian matrix $\mathcal{L}$ and its characteristic valuation as a function of $p$, given the single-layer network Laplacians $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$, they proved that

[TABLE]

Thus indicating that the algebraic connectivity of the interconnected system follows two distinct regimes, one in which its value is independent of the structure of the two layers, and the other in which its upper bound is limited by the algebraic connectivity of the weighted superposition of the two layers whose Laplacian is given by $\frac{1}{2}\left(\mathcal{L}_{A}+\mathcal{L}_{B}\right)$. More importantly, the discontinuity in the first derivative of $\lambda_{2}$ indicates a radical change of the structural properties of the system at $p^{*}$. Such dramatic change is visible in the coordinates of characteristic valuation of the nodes of the two network layers. In the regime $p\leq p^{*}$, the components of the vector corresponding to network $A$ are

[TABLE]

This means that the two network layers are structurally disconnected and independent. For $p\geq p^{*}$, we have

[TABLE]

which means that the components of the vector corresponding to interdependent nodes of network $A$ and $B$ have the same sign, while nodes in the same layer have alternating signs. Thus in this second regime, the system connectivity is dominated by inter-layer connections, and the two network layers are structurally indistinguishable.

Further studies Radicchi (2014) on the spectral properties of the normalized Laplacian of random interconnected networks show richer behavior than those of isolated networks. Depending on the combination of basic structural properties as e.g. degree distribution, degree correlation, and strength difference between intra- and interlayer connections, topological and dynamical phase transitions associated with the spectrum of the normalized Laplacian can be either discontinuous or continuous, and different regimes can disappear or even coexist. The rigorous generalization to the previous findings to m-layers and its application in diffusion processes can be found in Darabi Sahneh et al. (2015); Van Mieghem (2016); Cozzo and Moreno (2016).

V.3.2 Abrupt percolation transitions

In 2009-2010 the analysis of the percolation properties of interdependent networks started Leicht and D’Souza (2009); Buldyrev et al. (2010), with the influential paper by Buldyrev et al. (2010) addressing dependencies between networks. Focusing on the specific case of two interconnected networks A and B, they proposed a simple model of cascading failure propagation, and they discover a sharp transition in the probability of existence of a mutual giant component. The mutual giant component is defined as the set of nodes in network A which belong to an infinite-sized cluster in network A and also have their corresponding nodes in network B belong to an infinite-sized cluster in network B (or vice-versa). Their main focus in on the probability of the existence of a large cascading failure in the mutual giant component. The proposed model was analytically studied in terms of avalanches in Baxter et al. (2012) pointing out the existence of a hybrid phase transition. Several studies have focused on avoiding catastrophic failures via optimal coupling between networks based on their connectivity patterns e.g. Brummitt et al. (2012); Reis et al. (2014); Turalska et al. (2019), showing that interconnectivity can be an asset.

Later the problem in Buldyrev et al. (2010) was revisited by Son et al. (2012) demonstrating that the original study can be simplified conceptually by deleting all references to cascades of failures. The model is formulated in terms an order parameter $S$ which accounts for the probability that a randomly chosen node is a member of the infinite percolating cluster, thus they can apply classical percolation techniques. Nevertheless, the authors point out that for interdependent networks that are correlated with each other, or that are spatially embedded, the transition is in general not first order, and that interdependence can make the transition even less sharp than found in ordinary percolation.

The apparent mismatch between the results of Buldyrev et al. (2010) and Son et al. (2012) has been recently solved in a unifying framework proposed by Lee et al. (2016). The idea is to formalize the problem in terms of a hybrid phase transition, akin to that first reported in Baxter et al. (2012), for the order parameter $m(z)$ which is the size of the giant mutually connected cluster per node ($z$ refers to the average degree of the network). As defined in Sec. II.4.5 the order parameter for a hybrid transition is characterized as:

[TABLE]

where $m_{0}$ and $r$ are constants and $\beta_{m}$ is the critical exponent of the order parameter. Recall that the exponent must obey $0<\beta<1$ for the hybrid transition to hold. The first distinction offered by the authors in the cascade failure model, is between clusters and avalanches. Clusters are mutual connected components. Avalanches consist of mutual connected components separated from the giant component as a consequence of a removed node that triggers subsequent cascade. While the avalanche dynamics exhibits a critical pattern for a critical average degree $z_{c}$, the finite mutual connected components at $z_{c}$ mostly consist of one or two nodes, which is in discord with the power-law behavior of the cluster size distribution at a transition point characteristic of the conventional second-order percolation transition. From the proposed hybrid phase transition scheme, the authors classify avalanches in the critical region as finite and infinite avalanches. When an infinite avalanche occurs, the mutual giant component collapses, and the system falls into an absorbing state.

Apart from the previous abrupt transitions found, several authors have observed a sharp contrast to ordinary percolation between multilayer networks and traditional (monolayer) networks. The principal difference is rooted in the definition of the mutual component whereby a node is a member of the mutual component only if all the nodes that are interdependent with this node also belong to the mutual component. In monolayer networks, nodes of high degree belong to the giant connected component (percolation cluster) with higher probability, however, in multilayer networks the inclusion of nodes of high “superdegree” (that indicates the number of interdependent nodes in other layers of each individual node) into the mutual component is more difficult. The nontrivial point is that when the “superdegree” of the nodes in the layers are distributed heterogeneously among the layers, each given degree can enter the mutual component not smoothly but through a discontinuous transition Bianconi and Dorogovtsev (2014). It has been shown that for a large class of multilayer networks, the mutually connected component does not depend on the global topology. The mutual component problem for these networks can be solved by reduction to that for multiplex networks Bianconi et al. (2015). In multiplex networks, Cellai et al. (2013) presented a general framework for studying the emergence of a mutual giant component in multiplex networks with overlap. They show the presence of a critical value of edge overlap in a two-layer multiplex network that can change the order of the phase transition from hybrid first order to second order. The model is a directed version of percolation of interdependent multiplex networks with overlap which admits an epidemic spreading interpretation. This generalized percolation model reduces the to the model of interdependent multiplex networks in absence of overlap but in presence of overlap it is distinct. In Cellai et al. (2016) the authors provide a message passing theory for characterizing the percolation transition in multiplex networks with link overlap and an arbitrary number of layers, showing that the percolation transition has a discontinuous phase transition, with the exception of the trivial case in which all the layers completely overlap.

Recently, the authors in Radicchi and Bianconi (2017) propose a percolation model where the condition that makes a node functional is that the node is functioning in at least two of the layers of the network, redefining then the concept of mutual giant component in an M-layer multiplex setting. According to this model, interdependencies make the system more fragile than it would be by considering each layer independently. This fact is consistent with the traditional model used to study percolation in multiplex networks, and it is apparent from the fact that the percolation transition is abrupt for any number of network layers considered in the inter-dependent model.

The nuanced contrasts between the specific results reviewed here on interdependent networks show that the landscape is broad and that interconnectivity between networks can be a liability or an asset depending on the context.

V.3.3 Abrupt synchronization transitions

Another example of explosive phenomena in multilayer networks is that of explosive synchronization described in Sec.III. Several works have found abrupt synchronization transitions in this framework. We review them here.

In Zhang et al. (2015) the authors address the adaptive scheme described in Sec.III.7 in scale-free networks as layers of a multiplex network. In particular, the authors study a multiplex system composed of two networks (A and B) in which each of the $N$ nodes of a layer is associated (coupled) to another node of the other network. The evolution equations for the two layers read:

[TABLE]

where we denote with the same index $i$ each pair of nodes, one belonging to A and the other to B, that are coupled one-to-one. Matrices ${\bf A}$ and ${\bf B}$ correspond to the adjacency matrices of layer A and B respectively.

Given Eqs. (100)-(101), the nodes are split into the same two groups I and II as in the case of a single network. In this case, if a node $i$ in layer A is of type I (II) the corresponding node $i$ in layer B will be also of type I (II). Those $(1-f)N$ pairs of nodes of type II are not dynamically coupled, $\alpha_{i}^{A}=\alpha_{i}^{B}=1$, and they are only affected by the dynamics of their neighbors in their corresponding layer. On the other hand, if the pair of nodes are of type I they interact directly in the following way: $\alpha_{i}^{A}=r_{i}^{B}(t)$ and $\alpha_{i}^{B}=r_{i}^{A}(t)$. Therefore, the $fN$ nodes of type I in one layer are affected by the local degree of synchronization among the neighbors of node $i$ in the other network layer.

Figure 43 reports the synchronization diagrams for $f=1$ in three different ER-ER multiplex networks [panels (a), (b) and (c)], while panel (d) shows the transition obtained when the multiplex is composed by an ER and a SF layer with the same average degree $\langle k\rangle=12$. Interestingly, in all the four cases the explosive transition happens almost simultaneously for both layers.

A similar approach is developed in Huo et al. (2016), but they also impose a new phase to the Eqs. (100)-(101) in the multiplex defined as $\theta=p\theta^{(A)}+(1-p)\theta^{(B)}$. In the particular case of a power-law degree distributed network (A), and a Watts-Strogatz network (B), the weighted version of the phase $\theta$, undergoes a first order phase transition as a function of $p$. Yet, recently, Khanra et al. Khanra et al. (2018) also investigated a coupled phase-frustrated Sakaguchi-Kuramoto dynamics on a bilayer multiplex network. Essentially, the authors investigate Eqs. (100)-(101) in a framework where $\theta^{X}_{j}-\theta^{X}_{i}\rightarrow\theta^{X}_{j}-\theta^{X}_{i}-\alpha$, for every layer $X$. The results show that the frustrated system can undergo a sharp transition when coupled in a multilayer network for a wide region of the parameter space, moreover the critical points of the forward and backward transition can be computed analytically.

Very recently Nicolaou et al. Nicolau et al. (2019), inspired by the Kuramoto dynamics, introduced a new nonlinear model, called Janus oscillators, in which each unit $i$ splits into two oscillators. These oscillators have independent phases, $\theta_{i}^{A}$ and $\theta_{i}^{B}$, and different natural frequencies, $\omega_{i}^{A}$ and $\omega_{i}^{B}$, that are related in the following way: $\omega_{i}^{A}=\Omega+\omega/2$ and $\omega_{i}^{B}=\Omega-\omega/2$ (where as usual $\Omega$ is the average frequency of oscillators). These two faces of the same unit $i$ can be viewed as the two representations of each node $i$ in the layers of a duplex. However, the duplex induced by the coexistence of two oscillators in each dynamical unit $i$ takes a particular form, as revealed from the dynamical equations at work:

[TABLE]

From these equations we observe that the two oscillators in dynamical unit $i$ interact with each other through coupling $\lambda^{\prime}$ while oscillator $i$ in layer $A$ ($B$) interacts with those neighbors in layer $B$ ($A$) as dictated by the adjacency matrix ${\bf A}$. Thus, this duplex takes a bipartite form in the sense that the interactions of oscillators in layer $A$ ($B$) are always with those in layer $B$ ($A$). This structure automatically induces a frequency gap $\omega$ between each pair of connected oscillators. The frequency gap induced by the Janus model produces a variety of interesting behaviors, even in the case of simple network structure. Among these behaviors, the authors found ES by varying $\lambda$, thereby tuning from the asynchronous state to one of multiple synchronized states that can range from phase-locking to a variety of partially synchronized states.

Another interesting phenomena arises when coupling two different dynamics (one on each layer) interacting via a duplex (two-layer multiplex). Critical phenomena associated with this kind of interaction where first analyzed by Granell et al. (2013, 2014) in the context of epidemic spreading. However, the observation of an abrupt transition for interacting dynamical processes was found again in the context of synchronization. In Nicosia et al. (2017) the authors couple two dynamical processes, one for each layer of a duplex with their respective topologies. Assuming that the dynamics of the entire system are governed by the following equations:

[TABLE]

where ${\bf x}$ and ${\bf y}$ denote the vectors of states of the two dynamical processes, while the topologies of the two layers are encoded in the adjacency matrices $A$ and $B$ respectively.

The dynamical evolution of the two network processes are governed respectively by the functions $F_{\omega}$ and $G_{\chi}$, which depend on the sets of parameters $\omega$ and $\chi$, so that the state $x_{i}$ ($y_{i}$) of node $i$ at layer $A$ ($B$) is a function of the state ${\bf x}$ (${\bf y}$) and of the topology $A$ ($B$) of the first (second) layer. The key ingredient that connects the two dynamical processes is provided by the nature of the correspondence between layers. In fact, the parameter $\omega_{i}$ in function $F_{\omega_{i}}$ at layer A is itself a function of time which depends on the dynamical state $y_{i}$ of node $i$ at layer B, while the parameter $\chi_{i}$ at layer B depends on the state $x_{i}$ of node $i$ at layer A. Namely, we have:

[TABLE]

where $f$ and $g$ are two assigned functions. The specific problem analyzed consists of the coupling between a Kuramoto model in one layer, and a biased random walk on the network in the other layer:

[TABLE]

where $\omega_{i}$ corresponds to the natural frequency of the oscillator $i$ and $\lambda$ is the coupling strength of the Kuramoto model, and

[TABLE]

where $\pi_{ij}$ is the transition probability from node $j$ to node $i$, $\tau_{y}$ is the time scale of the random walker dynamics, and we have assumed that the random walk is biased on a node property $\chi_{i}$, with a tuneable bias exponent $\alpha$. The results of the model are depicted in Fig. 44, where the explosive synchronization phenomenon is revealed for the Kuramoto order parameter $r$. Very recently, the general formalism introduced in Eq. (104), has been used to study a multiplex in which a SIS epidemic model and the Kuramoto dynamics coexist Soriano-Paños (2019). The authors have found that both epidemic and synchronization transitions become discontinuous when the two dynamics interact in a cooperative way. This happens when synchronization in the first layer is promoted by the spreading activity in the second one and, simultaneously, the contagion rate in the second layer is fostered by the dynamical coherence in the first.

Finally, during the revision of this manuscript, a relevant result concerning the emergence of different types of synchronization transitions in multilayer networks has appeared Danziger et al. (2019). In this work, the authors analyze interdependent and competitive interactions between dynamical processes on multilayer networks, including synchronization dynamics. They define a function that accounts for the interaction between layers of the multiplex, in particular the influence that each node $i$ in layer A has on its replica node $i$ in layer B for a two-layered multiplex. The function relies on a local order parameter, $z_{i}$, on one layer, specified in Danziger et al. (2019) to be the local level of synchronization on that layer. This local order parameter is the input signal for the corresponding replica node in a different layer, for example: $F_{i}^{A\rightarrow B}(t)=|z_{i}^{A}|(t)$ specifies how much a node $i$ in layer A affects a node $i$ in layer B. For competitive interactions (i.e., those that discourage synchronization), the proposed form of the function is $F_{i}^{A\rightarrow B}(t)=|1-z_{i}^{A}|(t)$. Using this last specification, the global synchronization level undergoes a discontinuous backward transition as the coupling strength is decreased, reminiscent of an explosive process resulting from the frustration imposed on the replica nodes by that choice for the interaction function between layers.

VI Conclusions and Open Challenges

The paradigm of explosive phenomena in complex networks offers many novel advances, from the basic theoretical understanding of phase transitions to the ability to model a broad array of new phenomena on networks. Despite the large body of work on explosive phenomena and efforts to develop a deeper understanding, many abrupt transitions are still not well understood. This review offers an attempt to organize the extensive existing literature into thematic categories and identify fruitful directions for future research.

Since the start of vigorous activity in the field in 2009, there has been a chain of fascinating discoveries about the different exotic behaviors and universality classes of explosive phenomena. More specifically, we now understand a myriad of mechanisms that lead to EP and ES and have discovered many processes and models that implement these mechanisms. Through EP in particular we see the massive impact that a repeated, small intervention can achieve. This also provides new opportunities for modeling network failure modes in technological and economic systems, as well as interventions and consequences, as discussed briefly in Helbing (2013).

The phenomenology behind EP and ES share similar ingredients. In particular, both EP and ES are based on the implementation of some microscopic rules that avoid the formation of a macroscopic (connected or synchronized) component. In EP the rules are typically explicit and lead to the growth of a powder-keg distribution for the cluster-sizes at the percolation transition point. In ES the rules are implicit and result from the interplay of the natural frequency distribution and the network connectivity patterns. More understanding of the sub-critical synchronized cluster evolution for ES on different network topologies may provide an explicit framework to connect EP and ES.

Beyond the tour de force needed for establishing a mathematically grounded connection between EP and ES, there are still several challenges that remain unexplored for each phenomenon independently. For example, in EP, there is the need to quantify the tradeoff between the extent of global versus local information required. In ES, there is the question about how local order develops, e.g. do chimera states exist? EP has opened new opportunities for modeling networked systems, such as modular networks, optimal immunization strategies, network dismantling, and information spreading. Major outstanding challenges in the subject of EP are developing tailored control interventions for such systems. Major challenges in the subject of ES are that of relating the spectral properties of the graph with the outcome of the synchronization transition, or that of finding real-world examples where ES could be a plausible mechanism for switching on and off synchronization. ES is proving to be a useful paradigm for understanding neural disorders, which is extremely promising.

Many models of EP and ES have been designed to produce anomalous features via a built-in cluster growth suppression mechanism. Identifying additional natural, physical and socio-technological systems, where the suppression of large cluster mergers is an intrinsic feature of the dynamics is a compelling research avenue. The existence of such features in real-world systems may help to better understand and control these systems.

The array of behaviors, underlying mechanisms, and processes that implement explosive transitions on complex networks are still being discovered. And the explosive transitions surveyed in this review lead to many new research questions. We refer the reader to a number of recent reviews with extensive additional details and perspectives Bastas et al. (2014); Araújo et al. (2014); Saberi (2015); Boccaletti et al. (2016).

In summary, abrupt transitions may underlie many important phenomena and functions of systems found in nature, and it is of great importance to be able to describe and understand the mechanisms through which they occur. The explosive phenomena reported so far are opening up new perspectives on the types of transitions we can observe in networked systems and we envisage that they will catalyze increased future research activity on dynamical processes unfolding on top of complex networks.

Acknowledgements.

We gratefully acknowledge support from funding agencies and the reviewers of this manuscript. In particular RMD acknowledges support from the U.S. Army Research Office MURI Award No. W911NF-13-1-0340 and Cooperative Agreement No. W911NF-09-2-0053, and DARPA Award No. W911NF-17-1-0077. JGG acknowledges support from MINECO through project No. FIS2017-87519-P, and from the Departamento de Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo (FENOL group E-19). AA acknowledges support by Ministerio de Economía y Competitividad (grant FIS2015-71582-C2-1, PGC2018-094754-B-C21), Generalitat de Catalunya (grant 2017SGR-896), and Universitat Rovira i Virgili (grant 2017PFR-URV-B2-41), ICREA Academia and the James S. McDonnell Foundation (grant # 220020325). We benefited from critical feedback from G. Bianconi, B. Ziff, B. Khang, S. Gómez, C. Granell, M. Timme, Y. Fender, and L. Fender.