# Unified framework for hybrid percolation transitions based on microscopic
dynamics

Hoyun Choi  CTP and Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea  Y.S. Cho  Department of Physics, Jeonbuk
National University, Jeonju 54896, Korea  Raissa D’Souza  Department of
Computer Science and Department of Mechanical and Aerospace Engineering,
University of California, Davis, California 95616, USA  Santa Fe Institute,
Santa Fe, New Maxico 87501, USA  János Kertész  Department of Network and Data
Science, Central European University, Quellen strasse 51, 1100 Vienna, Austria
B. Kahng<EMAIL_ADDRESS>CCSS, KI for Grid Modernization, Korea
Institute of Energy Technology, Naju, Jeonnam 58217, Korea

###### Abstract

A hybrid percolation transition (HPT) exhibits both discontinuity of the order
parameter and critical behavior at the transition point. Such dynamic
transitions can occur in two ways: by cluster pruning with suppression of loop
formation of cut links or by cluster merging with suppression of the creation
of large clusters. While the microscopic mechanism of the former is understood
in detail, a similar framework is missing for the latter. By studying two
distinct cluster merging models, we uncover the universal mechanism of the
features of HPT-s at a microscopic level. We find that these features occur in
three steps: (i) medium-sized clusters accumulate due to the suppression rule
hindering the growth of large clusters, (ii) those medium size clusters
eventually merge and a giant cluster increases rapidly, and (iii) the
suppression effect becomes obsolete and the kinetics is governed by the Erdős-
Rényi type of dynamics. We show that during the second and third period, the
growth of the largest component must proceed in the form of a Devil’s
staircase. We characterize the critical behavior by two sets of exponents
associated with the order parameter and cluster size distribution, which are
related to each other by a scaling relation. Extensive numerical simulations
are carried out to support the theory where a specific method is applied for
finite-size scaling analysis to enable handling the large fluctuations of the
transition point. Our results provide a unified theoretical framework for the
HPT.

##  I Introduction

Hybrid phase transitions, sometimes referred to as mixed-order transitions,
exhibit features of both first- and second-order transitions at the same
transition point Chaikin and Lubensky (1995); Stanley (1971). Recently, hybrid
phase transitions have been observed in various models and systems, including
models of wetting Blossey and Indekeu (1995) and DNA denaturation Kafri  _et
al._ (2000), spin models with long-range interactions in one dimension Bar and
Mukamel (2014); Thouless (1969); Dyson (1971); Aizenman  _et al._ (1988), and
short-range competing interactions, such as the Ashkin-Teller model, in scale-
free networks at critical endpoints Jang  _et al._ (2015), glass and jamming
models Toninelli  _et al._ (2006), epidemic contagion models on complex
networks Janssen and Stenull (2016); Choi  _et al._ (2017a); Hasegawa and
Nemoto (2018), synchronization Pazó (2005); Gómez-Gardeñes  _et al._ (2011);
Coutinho  _et al._ (2013); Park and Kahng (2019); Basnarkov and Urumov (2008);
Park  _et al._ (2020), and percolation Panagiotou  _et al._ (2011); Buldyrev
_et al._ (2010); Baxter  _et al._ (2012); Dorogovtsev  _et al._ (2006); Choi
_et al._ (2021); Zhou  _et al._ (2014a); Schwarz  _et al._ (2006); Son  _et
al._ (2012); Choi  _et al._ (2017b) models. In addition, a hybrid phase
transition was experimentally demonstrated in a colloidal crystal Alert  _et
al._ (2017). However, the understanding of the microscopic mechanism
underlying this type of phase transition is far from being complete Araújo
_et al._ (2014); Dorogovtsev  _et al._ (2008); Lee  _et al._ (2018).

In percolation on a network, a geometric phase transition, the control
parameter $t=L/N$ is the density of occupied (active) links with $L$ being the
number of occupied links, and $N$ the number of nodes. The order parameter $m$
is the probability that a node belongs to a giant cluster as $N\to\infty$
Erdős and Rényi (1960); Achlioptas  _et al._ (2009); Bollobás (2001). In
Bernoulli percolation Stauffer and Aharony (1992), where the occupation
probability of the links is independent, there is usually (e.g., on regular
lattices in dimensions higher than one or on a complete graph) a second-order
or continuous transition from the non-percolating to the percolating phase
characterized by the emergence of a giant component. The transition is
accompanied by critical phenomena, like the power law behavior of
thermodynamic quantities near the critical point, which are characterized by
critical exponents, e.g., $\beta,\gamma,\nu$, for the order parameter, the
susceptibility, and the correlation length, respectively Stauffer and Aharony
(1992); Christensen and Moloney (2005). The percolation transition can also be
described in terms of the probability distribution;
$p_{s}=s^{1-\tau}\hat{f}(s/s^{*})$ with $s^{*}\sim|t-t_{c}|^{-1/\sigma}$,
where $p_{s}\equiv sn_{s}$ is the probability that a node belongs to a cluster
of size $s$, $n_{s}$ is the normalized number of clusters of size $s$, and
$\hat{f}$ is a scaling function. Of all the exponents, only two are
independent.

Here, we focus on hybrid percolation transitions (HPTs) on complex networks.
In an HPT, the order parameter is discontinuous at $t_{c}^{-}$ when the
transition point $t_{c}$ is approached from below; and it exhibits critical
behavior on the other side as $t\to t_{c}^{+}$ from above. This behavior is
expressed as follows:

$m(t)=\left\\{\begin{array}[]{lr}0&~{}{\rm for}~{}~{}t<t_{c},\\\
m_{c}+r(t-t_{c})^{\beta_{m}}&~{}{\rm for}~{}~{}t\geq t_{c},\end{array}\right.$
(1)

where $m_{c}=m(t_{c})$; $r$ is a $t$-independent constant; and $\beta_{m}$ is
a critical exponent associated with the order parameter.

Bernoulli percolation can be considered via two possible methods: Either the
links are occupied with a probability corresponding to $t$ (static picture),
or the links are occupied one by one until their density reaches $t$ (kinetic
picture D’Souza and Nagler (2015); Ben-Naim and Krapivsky (2005)) - and these
two are equivalent in the thermodynamic limit. Here we are dealing with
dynamic models with correlated occupations, so we have to apply the kinetic
representation, which has the additional advantage that we can trace the
evolution of the characteristic quantities.

HPTs can occur both in cluster pruning ($t\to t_{c}^{+}$) and cluster merging
($t\to t_{c}^{-}$) processes. For cluster pruning ($t\to t_{c}^{+}$), it is
established that the discontinuity of the order parameter is related to a
dynamically evolving metastable state while the critical behavior is induced
by cascading failures with the underlying universal mechanism identified as a
critical branching process Lee  _et al._ (2017); Zhou  _et al._ (2014b). To
characterize the critical behaviors, two sets of critical exponents were
introduced Lee  _et al._ (2016a): one is associated with the order parameter
(i.e., the giant cluster) and the other is associated with the power-law
behavior of the avalanche size distribution; moreover, a scaling law relates
these two sets of critical exponents.

For the case of the cluster-merging process such a universal mechanism has yet
to be identified. A challenge is that the dynamic evolution of the system
starts from $t=0$, thus the critical behavior of HPT is affected by the
process leading to the jump of the order parameter. The following natural
questions arise: Is there a universal microscopic mechanism behind the HPT in
cluster merging processes? Can the critical behaviors of the order parameter
and the clusters be characterized by critical exponents like for the cluster
pruning processes? Is there a relationship between the critical behaviors of
the cluster pruning and the cluster merging processes? These are the questions
we address in this paper.

The paper is organized as follows. In Sec. II, we review the results on the
microscopic mechanism of the critical phenomena in the HPT due to cluster-
pruning processes. In Sec. III.1 and III.2 we reconsider the critical behavior
of the $r$-ER model, and introduce and consider the $m$-BFW model from a
common perspective; we discuss the universal scaling relation that unifies the
critical behaviors in the cluster-merging and cluster-pruning processes. In
Sec. IV, we summarize the results and discuss their implications.

##  II HPT in cluster-pruning process

In this section we present the established results for the HPT in the cluster
pruning process Lee  _et al._ (2016a), for instance, the critical behavior of
$k$-core percolation with $k=3$ on the Erdős-Rényi (ER) graph Erdős and Rényi
(1960). In this case, the kinetics starts from a supercritical state at $t\gg
t_{c}$. The $k$-core Chalupa  _et al._ (1979); Baxter  _et al._ (2015); Lee
_et al._ (2016b) is a subgraph of the network in which each node has at least
$k$ degrees. To obtain an initial $k$-core subgraph, an ER network is
generated at $t\gg t_{c}$, and then all nodes with a degree of less than $k$
are iteratively removed along with their links. From the resulting stable
$k$-core subgraph, a randomly selected node is removed along with its links.
This removal may lead to other nodes getting knocked out of the $k$-core in an
avalanche of activity, reducing the size of the $k$-core. The amount by which
the size is reduced corresponds to avalanche size. As this process is
repeated, the order parameter decreases gradually and exhibits critical
behavior following Eq. (1). At $t_{c}^{+}$, the order parameter suddenly drops
from a finite value to zero.

The features of this critical behavior are as follows. (i) The critical
exponent can be derived analytically as $\beta_{m}=1/2$ Lee  _et al._ (2016a),
which is universal regardless of the different models, such as the cascade
failure model of interdependent networks Buldyrev  _et al._ (2010). (ii) The
critical behavior for $t>t_{c}$ is characterized by two sets of critical
exponents: one set $\\{\beta_{m}$, $\gamma_{m}$, $\bar{\nu}_{m}\\}$ is
associated with the order parameter, and the other set
$\\{\tau_{a},\sigma_{a},\beta_{a},\gamma_{a},\bar{\nu}_{a}\\}$ is associated
with the avalanche size distribution. For instance, the fluctuations of the
order parameter over different configurations are expressed by the
susceptibility, which is defined as $\chi_{m}\equiv N(\langle
m^{2}\rangle-\langle m\rangle^{2})\sim(t-t_{c})^{-\gamma_{m}}$. However, the
mean avalanche size is expressed as
$\chi_{a}\equiv\sum_{s}^{\prime}sp_{s}(t)/\sum_{s}^{\prime}p_{s}(t)\sim(t-t_{c})^{-\gamma_{a}}$,
where the prime denotes the summation over finite avalanche sizes, and
$p_{s}(t)$ denotes the avalanche size distribution at $t$. As
$\gamma_{m}\approx 1.0$ and $\gamma_{a}=1/2$, the two $\gamma$-s are different
(though respectively universal for the different cluster pruning models on the
ER model Lee  _et al._ (2016a)). The exponents $\tau_{a}$, $\sigma_{a}$, and
$\bar{\nu}_{a}$ characterize the avalanche size distribution: $n_{s}\sim
s_{a}^{-\tau_{a}}f(s_{a}/s_{a}^{*})$, where $s_{a}$ denotes the avalanche
size, $f$ is a scaling function, and $s_{a}^{*}$ is the characteristic
avalanche size, which behaves as $\sim(t-t_{c})^{-1/\sigma_{a}}$ for
$N\to\infty$. Its finite size behavior at the transition point is given by
$s_{a}^{*}\sim N^{1/\sigma_{a}\bar{\nu}_{a}}$. (iii) The critical exponents of
the individual set satisfy the scaling relations; however, the two sets are
not independent, but they are related through a conservation relation Lee  _et
al._ (2016a):

$m(t)+\int_{t}^{t_{0}}\sum_{s}^{\prime}sp_{s}(t_{1})dt_{1}=1,$ (2)

leading to

$\gamma_{a}=1-\beta_{m}.$ (3)

This scaling relation is universal for the cascading failure models. Later we
will show that the relation $\eqref{eq:gamma_a_beta_m}$ is also valid for the
HPT in the cluster-merging process.

##  III HPT in cluster-merging process

###  III.1 Modified $r$-ER model

$g$ |  $\tau_{s}$ |  $\sigma_{s}$ |  $\beta_{s}$ |  $\bar{\nu}_{s}$ |  $\gamma_{s}$ |  $\gamma_{s}^{\prime}$
---|---|---|---|---|---|---
0.2 |  $2.08\pm 0.04$ |  $0.99\pm 0.05$ |  $0.09\pm 0.05$ |  $1.10\pm 0.10$ |  $0.91\pm 0.05$ |  $1.03\pm 0.02$
0.5 |  $2.18\pm 0.04$ |  $0.96\pm 0.05$ |  $0.19\pm 0.05$ |  $1.24\pm 0.11$ |  $0.83\pm 0.05$ |  $1.10\pm 0.03$
0.8 |  $2.25\pm 0.04$ |  $0.89\pm 0.05$ |  $0.29\pm 0.07$ |  $1.42\pm 0.13$ |  $0.81\pm 0.05$ |  $1.15\pm 0.05$
$g$ |  $\beta_{m}$ |  $\bar{\nu}_{m}$ |  $\gamma_{m}$ |  $\gamma_{m}^{\prime}$ |  $\zeta$ |  $\zeta^{\prime}$
0.2 |  $0.09\pm 0.05$ |  $1.05\pm 0.10$ |  $0.90\pm 0.10$ |  $2.06\pm 0.05$ |  $1.04\pm 0.05$ |  $1.89\pm 0.10$
0.5 |  $0.21\pm 0.05$ |  $1.25\pm 0.15$ |  $0.83\pm 0.10$ |  $2.30\pm 0.10$ |  $1.02\pm 0.05$ |  $1.76\pm 0.13$
0.8 |  $0.32\pm 0.05$ |  $1.40\pm 0.15$ |  $0.76\pm 0.10$ |  $2.55\pm 0.10$ |  $1.00\pm 0.07$ |  $1.65\pm 0.13$
Table 1:  List of numerical values of the exponents for HPT in the modified
$r$-ER model. The numerical values of all critical exponents except $\zeta$
and $\zeta^{\prime}$ are newly obtained. They are consistent within the error
bars with those presented in Ref. Cho  _et al._ (2016). The exponents values
of $\zeta$ and $\zeta^{\prime}$ are adopted from Ref. Park  _et al._ (2019).
The definitions of each critical exponent and the scaling relations among them
are listed in Appendix A.

We next present results for the HPT due to cluster merging for two distinct
models and establish the underlying universal mechanism. When HPT occurs as a
consequence of a cluster merging process the control parameter moves in the
opposite direction as compared to cluster pruning processes from small to
large values of $t$: the order parameter first jumps, and then the critical
behavior is observed. Thus, the critical behavior is affected by the dynamics
that occur leading to the jump of the order parameter as $t\to t_{c}^{-}$. The
first model we study is the restricted ER model (denoted as the $r$-ER model).
This model was originally proposed to study a discontinuous percolation
transition Panagiotou  _et al._ (2011); here we consider a slightly modified
version, the modified $r$-ER model, which was shown to exhibit an HPT Cho  _et
al._ (2016) (see Appendix B).

In this model, the cluster coalescence dynamics begin with $N$ isolated nodes.
At each time step, the clusters are ranked by size and partitioned into one
set of small and one set of large clusters, denoted as $A$ and $B$,
respectively. In the original $r$-ER model, the $g$ fraction of nodes
contained in the smallest clusters is assigned to set $A$, where $g$ ($0<g<1$)
is a model parameter. The remaining fraction of nodes is assigned to set $B$.
Hence one cluster can have nodes belonging to both sets (see Fig. A1). We make
the modification that when a cluster has nodes that span both sets $A$ and
$B$, all nodes in the cluster are regarded as elements of set $A$. From there
we proceed as with the original $r$-ER model. Two nodes are randomly selected
from different clusters for cluster coalescence: one is from the entire system
and the other is only from set $A$. Thus, two nodes in set $B$ cannot be
linked and this restriction rule suppresses the growth of large clusters.
Accordingly, large clusters are rarely generated, whereas medium-sized
clusters are generated in abundance. As $t$ approaches $t_{c}$, medium-sized
clusters are more likely to be merged, and the large cluster size rapidly
increases. In the limit $N\to\infty$, the growth rate of the giant cluster
with respect to $t$ becomes infinity, and the change of the order parameter
becomes discontinuous.

In finite systems, the cluster-merging process is characterized by three time
regimes: $[0,t_{a}]$, $[t_{a},t_{g}]$, and $[t_{g},t_{c}]$ Park  _et al._
(2019) (see Appendix B). In the early regime $t\in[0,t_{a}]$, due to the
suppression effect, a bump is formed in the cluster-size distribution and its
size increases first with $t$ corresponding to a large number of medium-size
clusters. Here $t_{a}$ is defined as the time at which the bump size reaches
its maximum. The bump is consistent with a “powder keg” as first discussed by
Friedman and Landsberg Friedman and Landsberg (2009). In the intermediate
regime $[t_{a},t_{g}]$, the bump size shrinks as medium-sized clusters merge
with the giant cluster, and the size of the giant cluster increases rapidly.

At $t_{g}$, partition $B$ is fully occupied by the giant cluster alone; that
is, $m(t_{g})=1-g$. Hence, for $t\geq t_{g}$, the giant component must contain
nodes in set $A$, thus the partition into two categories loses its meaning. In
the regime $[t_{g},t_{c}]$, the cluster-merging kinetics proceeds according to
the ER rule Cho  _et al._ (2010); however, the initial configuration of the ER
kinetics is not the standard initial condition of isolated nodes, but the
specific cluster configuration at $t_{g}$. The size distribution of finite
clusters exhibits a power-law decay in a small cluster-size region with
exponent $\tau_{s}$, which is not the value of the ER (or mean-field) class
$\tau_{s}^{\rm{ER}}=5/2$ but depends on the model parameter $g$. When the
system reaches $t_{c}$ the bump is completely eliminated, and finite clusters
have a size distribution following a Stauffer-type cluster scaling form with
exponent $\tau_{s}$, i.e., $n_{s}\sim s^{-\tau_{s}}\exp(-s/s^{*})$, where
$s^{*}$ is a characteristic cluster size due to the finite-size effect:

$s^{*}\sim N^{1/\sigma_{s}\bar{\nu}_{s}},$ (4)

with system size $N$. The cluster scaling formula is valid also for
$t_{g}<t<t_{c}$ with

$s^{*}\sim(t_{c}-t)^{-1/\sigma_{s}^{\prime}},$ (5)

as long as $s^{*}$ is smaller than the value given by (4). The critical
behavior of the clusters can be described using exponents
$\\{\tau_{s},\sigma_{s},\bar{\nu}_{s}\\}$. Since the merging of finite
clusters continues to follow the ER mechanism also for $t>t_{c}$, the critical
behavior in this regime is characterized by two sets of critical exponents:
one set $\\{\beta_{m}$, $\gamma_{m}$, $\bar{\nu}_{m}\\}$ is associated with
the order parameter, and the other set
$\\{\tau_{s},\sigma_{s},\beta_{s},\gamma_{s},\bar{\nu}_{s}\\}$ is associated
with the size distribution of finite clusters.

Figure 1:  Measurement of the critical exponents of the modified $r$-ER model
for $g=0.2$. (a) Plot of $n_{s}(m_{c})$ vs $s$. $n_{s}(m_{c})$ is measured
just after $m$ exceeds $m_{c}$ for different system sizes
$N/10^{4}=2^{0},2^{2},2^{4},2^{6},2^{8},~{}{\rm and}~{}2^{10}$ from the left.
Here, $m_{c}$ of each $N$ is taken as a value immediately after the bump of
$n_{s}(m_{c})$ disappears. A dashed line with slope $-2.08$ is drawn to guide
the eye. (b) Scaling plot of $s^{\tau_{s}}n_{s}(m_{c})$ vs $s/s^{*}$ for
different sizes $N$ in (a), where $\tau_{s}=2.08$ and
$s^{*}=N^{1/\sigma_{s}\bar{\nu}_{s}}$ with $\sigma_{s}=0.99$ and
$\bar{\nu}_{s}=1.10$. (c) Plot of $t_{c}(N)-t_{c}(\infty)$ vs $N$ with error
bars for $t_{c}(N)>t_{c}(\infty)$, where $t_{c}(N)$ is the configuration
average of $t$ just after $m(t)$ exceeds $m_{c}(N)$. The slope of the solid
line is $-1/1.05$, which supports that $\bar{\nu}_{m}\approx 1.05$ by the
relation $t_{c}(N)-t_{c}(\infty)\sim N^{-1/\bar{\nu}_{m}}$. (d) Plot of
$m_{c}(N)-m_{c}(\infty)$ vs $N$ with error bars. The slope of the solid line
is $-0.07$, which supports that $\beta_{m}/\overline{\nu}_{m}\approx 0.07$ by
the relation $m_{c}(N)-m_{c}(\infty)\sim N^{-\beta_{m}/\bar{\nu}_{m}}$. (e)
Plot of $\chi_{m}(t)$ vs $(t-t_{c}(\infty))$ for $t>t_{c}(\infty)$ using
$N/10^{4}=2^{0},2^{2},2^{4},2^{6},2^{8},~{}{\rm and}~{}2^{10}$ from the right,
where $\chi_{m}(t)=N(\langle m(t)^{2}\rangle-\langle m(t)\rangle^{2})$. The
slope of the dashed line is $-0.90$, which supports that $\gamma_{m}\approx
0.90$ by the relation $\chi_{m}(t)\sim(t-t_{c}(\infty))^{-\gamma_{m}}$. (f)
Plot of $N\chi_{m}(t)$ vs $(t_{c}(\infty)-t)$ for $t<t_{c}(\infty)$ using
$N/10^{4}=2^{0},2^{2},2^{4},2^{6},2^{8},~{}{\rm and}~{}2^{10}$ from the right.
The slope of the dashed line is $-2.06$, which supports that
$\gamma_{m}^{\prime}\approx 2.06$ by the relation
$N\chi_{m}(t)\sim(t_{c}(\infty)-t)^{-\gamma_{m}^{\prime}}$.

We note that even though the evolution proceeds under the ER rule in regime
$[t_{g},t_{c}]$, the critical exponent $\tau_{s}$ does not have the ER value.
Instead, it is determined by the early time process. As clusters are merged
and their rankings are updated, they may move back and forth between the two
sets. The inter-event time (denoted as $z_{i}$) is defined as the time
interval between two set-changing events of the cluster containing node $i$
Park  _et al._ (2019). Its distribution $P_{I}(z)$ comprising these inter-
event times over all nodes during the time interval $[0,t_{a}]$ exhibits a
power-law decay, $P_{I}(z)\sim z^{-\zeta}$. It was derived Park  _et al._
(2019) that the exponent $\zeta$ is related to the exponents $\tau_{s}$ and
$\sigma_{s}$ as $\zeta=4-(\tau_{s}+\sigma_{s})$. Because the ER dynamics are
dominant in $[0,t_{a}]$ and the ER values are $\tau_{s}=5/2$ and
$\sigma_{s}=1/2$, $\zeta=1$ is predicted. Indeed, numerical simulations
produce $\zeta=1$, which is insensitive to model parameter $g$.

The exponent of $P_{I}(z)$ during the time interval $[t_{a},t_{g}]$ is denoted
by $\zeta^{\prime}$. Since in this interval, $s^{*}$ can be considered as
constant, $\sigma_{s}^{\prime}=0$. Thus, the above relationship changes to
$\zeta^{\prime}=4-\tau_{s}$ Park  _et al._ (2019). The dependence of
$\zeta^{\prime}$ on $g$ indicates that the critical exponent $\tau_{s}$ also
depends on it.

We present numerical values of the two sets of critical exponents for
different values of $g$ in Table. 1. The critical exponents
{$\tau_{s},\sigma_{s},\gamma_{s},\gamma_{s}^{\prime},\beta_{m}$} for different
$g$ were directly measured from numerical simulations. Then, $\bar{\nu}_{s}$
and $\beta_{s}$ are calculated using the scaling relations
$\beta_{s}=(\tau_{s}-2)/\sigma_{s}$ and
$\bar{\nu}_{s}=(\tau_{s}-1)/\sigma_{s}$, respectively. When the system reaches
$t_{c}$ ($m$ reaches $m_{c}$), we check that $s^{\tau_{s}}n_{s}$ vs
$s/s*=sN^{-1/\sigma_{s}\bar{\nu}_{s}}$ for different system sizes are well
collapsed using $\sigma_{s}$ and $\bar{\nu}_{s}$ within the error range of
numerical values, irrespective of $g$ (See Fig. 1(a) and (b) for the result of
$g=0.2$.). $\bar{\nu}_{m},\gamma_{m}$, and $\gamma_{m}^{\prime}$ are obtained
from simulation data as shown in Fig. 1(c)–(f). The values of $\zeta$ and
$\zeta^{\prime}$ are taken from Park  _et al._ (2019).

Next, we check whether the numerical values of exponents in Table. 1 satisfy
the scaling relations. For $g=0.2$, the exponent of the order parameter
$\beta_{m}$ defined in the formula $m-m_{c}\sim(t-t_{c})^{\beta_{m}}$ is
directly measured as $\beta_{m}=0.09\pm 0.05$. This value is consistent with
$\beta_{s}=0.09\pm 0.05$ obtained using the relation
$(\tau_{s}-2)/\sigma_{s}$. Secondly, using $\tau_{s}$ and $\sigma_{s}$,
$\gamma_{s}$ is obtained using the scaling relation $(3-\tau_{s})/\sigma_{s}$
as $\gamma_{s}=0.93\pm 0.09$. This value is consistent with the value
$\gamma_{s}=0.91\pm 0.05$ directed measured in Cho  _et al._ (2016) using the
formula of the average size of finite clusters $\langle
s\rangle\sim(t-t_{c})^{-\gamma_{s}}$, where $\langle
s\rangle\equiv\sum^{\prime}_{s}s^{2}n_{s}(t)/\sum^{\prime}_{s}sn_{s}(t)$.
Therefore, the scaling relation (3) between the two types of critical
exponents is valid within the error range as $\gamma_{s}+\beta_{m}=1.00\pm
0.10$. Note that we used $\gamma_{s}$ instead of $\gamma_{a}$. Finally, the
two exponents related to the inter-event time distribution $\zeta=1.04\pm
0.05$ and $\zeta^{\prime}=1.89\pm 0.10$ satisfy the relations $\zeta=1$ and
$\zeta^{\prime}=4-\tau_{s}=1.92\pm 0.04$ within the error ranges. We also
check the scaling relations for $g=0.5$, finding that they are valid within
the error ranges. However, for $g=0.8$, $\gamma_{s}+\beta_{m}=1.13\pm 0.05$
slightly deviates from unity. This discrepancy is caused by the deviation of
the estimated value $\sigma_{s}=0.89\pm 0.05$ from the theoretical value
$\sigma_{s}=1$ Cho  _et al._ (2016). We think that the origin of the deviation
is the crossover from HPT to continuous PT as $g\rightarrow 1$. Finally, we
confirm that the two scaling relations $\bar{\nu}_{s}=2\beta_{s}+\gamma_{s}$
and $\bar{\nu}_{m}=2\beta_{m}+\gamma_{m}$ are valid within the error ranges
using the values in Table. 1. However, we found that
$\chi_{m}(t)\propto(t-t_{c})^{-\gamma_{m}}$ is not clearly shown for $g=0.8$
again because of the crossover from HPT to continuous PT as $g\rightarrow 1$
as mentioned above. In conclusion, we demonstrated that the scaling relations
in HPT in cluster-pruning processes still hold in HPT in cluster-merging
processes by using the modified $r$-ER model.

###  III.2 Modified BFW model

we next study a fundamentally different model of the HPT, which does not rely
on partitioning the system into sets, and identify the same features as
obtained for the $r$-ER model.

To test the universality of the picture presented in III.1 we study a modified
version of the model proposed by Bohman et al (referred to as the $m$-BFW
model) Bohman  _et al._ (2004); Chen and D’Souza (2011), which does not
contain any partitions in the system, instead the suppression effect arises in
a self-organized manner via link rejection and a dynamic cap on the maximum
size component allowed.

The $m$-BFW model is characterized by three parameters, denoted by $k$ the
cap, $u$ the number of attempted link additions, and $h$ the lower bound on
the fraction of attempted link additions that must be accepted. While $h$ is a
constant in the range$(0,1]$, $k$ and $u$ evolve in time. The model is defined
as follows:

Figure 2:  Plots of the cluster size distributions of $m$-BFW model at (a)
$t\ll t_{a}$, (b) $t=t_{a}$, (c) $t=t_{c}$, and (d) $t>t_{c}$. Inset:
Schematic plots of the order parameter $m(t)$ vs $t\equiv L/N$ with a point
($\star$) indicating the corresponding time. (e) Plot of the order parameter
with characteristic points $\\{t_{a},t_{b},t_{c}\\}$ in the evolution of a
single sample of the $m$-BFW model with a system size $N=1.0\times 10^{4}$ for
$h=0.2$.

  * (i)

First, the system is initialized with $N$ isolated nodes, and no links, $L=0$.
Initially, $k(0)=2$ and $u(0)=0$ are set.

  * (ii)

Two nodes are selected randomly from two different clusters of size $s_{1}$
and $s_{2}$ and the number of trials is updated $u\to u+1$. A link between the
selected nodes is formed and that link is added ($L\to L+1$) if

    * (ii-1)

$s\equiv s_{1}+s_{2}\leq k$

    * (ii-2)

If $s\equiv s_{1}+s_{2}>k$ but $L/u\leq g$, then $k\to s_{1}+s_{2}$ is set.

Otherwise no link is formed and $L,k$ remain unaltered.

  * (iii)

Return to step (ii) and repeat the loop until a single cluster of size $N$ is
left.

In the original BFW model, $h(k)$ depends on $k$ as $h(k)=1/2+(2k)^{-1/2}$.
Instead here for the $m$-BFW model, $h$ is constant. This difference yields an
HPT. For details, see Appendix D. In this model, $k$ controls the size of the
largest cluster for a given $L$. The number of trials to add links is counted
by $u$. Thus, $L/u$ indicates the acceptance rate of the link attachment for a
given $u$, which is lower bounded by $h$. For convenience, we use $t\equiv
L/N$ as time (corresponding to the so-called event-time). Note that for $h=1$,
the $m$-BFW model is reduced to the ER model.

Large clusters in the $m$-BFW model are subject to two competing factors:
suppression of growth by the model rule, and higher growth rate, proportional
to their size. Accordingly, the largest cluster grows in a boom-and-bust
manner resulting in a Devil’s staircase pattern. This process effectively
performs a role similar to the back-and-forth dynamics across the set boundary
in the $r$-ER model. In fact, the cluster evolution of the $m$-BFW model is
similar to that of the $r$-ER model. We find that there also exist three time
intervals in the $m$-BFW model, denoted by $[0,t_{a}]$, $[t_{a},t_{b}]$, and
$[t_{b},t_{c}]$. The characteristic times $t_{a}$, $t_{b}$, and $t_{c}$ depend
on $N$ as shown in Fig. 2. In each time interval, the clusters evolve in a
manner similar to that in the corresponding regime of the $r$-ER model. In the
early time regime $t\in[0,t_{a}]$, a bump containing a large number of medium-
sized clusters forms, and its size increases in the cluster-size distribution.
In the intermediate time regime $[t_{a},t_{b}]$, the bump shrinks and the
giant cluster grows rapidly. In the late time regime $[t_{b},t_{c}]$, the bump
is eroded and a pure power-law cluster size distribution is built up at
$t_{c}$. In the late regime, the cluster merging proceeds via the ER dynamics,
but with the cluster configuration $n_{s}(t_{b})$ as an initial condition.
Note that the characteristic time $t_{b}$ corresponds to $t_{g}$ of the $r$-ER
model. $t_{g}$ is determined explicitly as $m(t_{g})=1-g$ in the $r$-ER model;
however, $t_{b}$ is determined in a self-organized manner.

Figure 3:  Plot of the order parameter $m(t)$ vs $t\equiv L/N$ in the $m$-BFW
model for $h=1.0$ (dashed, the ER case), $0.7$ (dotted), and $0.5$ (solid).

Fig. 3 shows a plot of $m$ vs $t$ for different values of the model parameter
$h$. We find that the $m$-BFW model exhibits HPTs for $h<1$, but a continuous
transition for $h=1$.

Figure 4:  (a) Acceptance rate $L/u$ (left axis) and component size cap $k$
(right axis) vs time $t$ for $h=0.2$ in the $m$-BFW model. The acceptance rate
$L/u$ exhibits increasing and then decreasing behavior on each plateau of $k$
as $t$ increases. The size cap $k$ exhibits a Devil’s staircase pattern for a
single configuration. Inset: Zoomed-in plot of the main panel in a given
region. (b) Plot of the mean step width of the Devil’s staircase
$\langle\Delta t\rangle$ vs $t$ in the form of $t_{c}(N)-t$. The regimes
$[0,t_{a}]$, $[t_{a},t_{b}]$, and $[t_{b},t_{c}]$ are shaded in blue, green,
and red (from right to left), respectively. The data are obtained from
$10^{5}$ configurations with a system size of $1.024\times 10^{7}$.

In Fig. 4(a), we plot the acceptance rate $L/u$ and the growth limit $k$ vs
$t$. We show that $k$ increases in a Devil’s staircase pattern. For a given
staircase plateau, the acceptance rate $L/u$ begins to increase and then it
decreases. The first part is because $k$ is upgraded each time $L/u$ decreases
to $h$, thus links can be more readily attached and the acceptance rate $L/u$
increases. As the size of the largest cluster increases and reaches the limit
$k$, the number of successful link attachment attempts declines repeatedly,
which reduces $L/u$ back to $h$. During this process, $k$ remains constant,
therefore, $k$ increases following a Devil’s staircase pattern. We plot the
average width of the staircases $\langle\Delta t\rangle$ over different
configurations as a function of $t_{c}(N)-t$ in Fig. 4(b). $\langle\Delta
t\rangle$ decreases in a power-law manner for $t<t_{a}(N)$, decreases rapidly
in $[t_{a}(N),t_{b}(N)]$, and is almost constant in $[t_{b}(N),t_{c}(N)]$.
Thus, $t_{b}(N)$ can be regarded as a crossover point to the ER limit, and for
$t_{b}(N)<t<t_{c}(N)$, $\langle\Delta t\rangle\approx 1/N$. This means that
almost every trial of the link attachment was accepted, particularly for giant
clusters without any restrictions and the dynamics of the link attachment
follows that of the ER model.

Figure 5:  (a) Plot of the cluster-size distribution $n_{s}(m;N)$ vs cluster
size $s$ for the $m$-BFW model with $h=0.2$. Data are collected from each
configuration when the giant cluster size first became larger than $mN$.
$m_{c,1}(N)$ is chosen when the bump in the cluster size distribution
disappears completely, and $n_{s}(m_{c,1};N)\sim
s^{-\tau_{s}}\exp{(-s/s^{*})}$ with $\tau_{s}\approx 2.04$. $s^{*}$ is the
characteristic size due to the finite-size effect. The area
$\sum_{s}^{\prime}n_{s}(m;N)$ (colored in blue (gray)) is denoted by
$1-m_{c,2}(N)$, which indicates the fraction of nodes contained in all finite
clusters. Accordingly, the order parameter was obtained as $m_{c,2}(N)$.
Simulations were performed in systems of size $N=1.024\times 10^{7}$ with the
model parameter $h=0.2$. Data are averaged over $10^{5}$ configurations. (b)
Plot of the finite-size scaling analysis. The characteristic size $s^{*}\sim
N^{1/{\sigma_{s}\bar{\nu}_{s}}}$ with $\sigma_{s}\bar{\nu}_{s}\approx 1.20.$

Due to the complexity of the transition, finite size scaling should be studied
with special care. We determine a transition point $t_{c}(N)$ as follows: We
denote the sample index by $i$. At the first time when the largest cluster
size is equal to or larger than $m$ the size distribution of finite clusters
divided by $N$ for sample $i$ is $n_{s}^{(i)}(m,N)$. We take an average over
all samples:

$n_{s}(m;N)\equiv\langle n_{s}^{(i)}(m;N)\rangle_{i}.$ (6)

When $m$ is small, the size distribution $n_{s}(m;N)$ of finite clusters
contains a bump. However, there exists a characteristic value $m_{c,1}(N)$ at
which the bump disappears completely, and the distribution $n_{s}(m_{c,1};N)$
exhibits power-law behavior. Therefore, $n_{s}(m_{c,1};N)$ can be expressed as
$\sim s^{-\tau_{s}}\exp(-s/s^{*})$, as shown in Fig. 5(a). We directly
measured the exponent $\tau_{s}$ as $\tau_{s}=2.04\pm 0.02$ for $h=0.2$. The
characteristic cluster size $s^{*}$ behaves as $(t_{c}-t)^{-1/\sigma_{s}}$ for
$t<t_{c}$ and $\sim N^{1/\sigma_{s}\bar{\nu}_{s}}$ at $t_{c}$. We estimate
that $1/\sigma_{s}\bar{\nu}_{s}\approx 0.83$ ($\sigma_{s}\bar{\nu}_{s}\approx
1.20$), as shown in Fig. 5(b). We define $m_{c,2}(N)\equiv
1-\sum^{\prime}_{s}sn_{s}(m_{c,1};N)$, where the prime denotes summation over
finite clusters. Next, for each sample $i$, we find $t_{c}^{(i)}(N)$, where
the giant cluster size is increased to $m_{c}^{(i)}(N)\geq m_{c,2}$ for the
first time. Subsequently, $t_{c}(N)$ is determined as $\langle
t_{c}^{(i)}(N)\rangle_{i}$, and $m_{c}(N)$ is determined as $\langle
m_{c}^{(i)}(N)\rangle_{i}$. Based on $t_{c}(N)$, we determine $t_{c}(\infty)$
using the finite-size scaling theory, and the plot of $t_{c}(\infty)-t_{c}(N)$
vs $N$ is assumed to follow a power law as $\sim N^{-1/\bar{\nu}_{m}}$ by
choosing $t_{c}(\infty)$ (see Fig. 6(a)). We estimate $t_{c}(\infty)=0.9953$
and $1/\bar{\nu}_{m}\approx 0.85$ ($\bar{\nu}_{m}\approx 1.18$) for $h=0.2$.
We also determine $m_{c}(\infty)$ using finite-size scaling:
$m_{c}(N)-m_{c}(\infty)\sim N^{-\beta_{m}/\bar{\nu}_{m}}$ by choosing the
proper $m_{c}(\infty)$ showing a power-law decay (see Fig. 6(b)). We obtain
$m_{c}(\infty)\approx 0.9555$ and $\beta_{m}/\bar{\nu}_{m}=0.05$ for $h=0.2$.
Therefore, $\beta_{m}\approx 0.062$ is estimated. We also directly measured
the exponent $\beta_{m}$ as $0.06\pm 0.01$ (see Fig. 8).

Figure 6:  Finite size scaling plots for the m-BFW model with $h=0.2$. (a)
$t_{c}(\infty)-t_{c}(N)$ vs $N$ with error bars for $t<t_{c}(\infty)$. It
behaves as $\sim N^{-1/\bar{\nu}_{m}}=N^{-0.85}$. (b) Plot of
$m_{c}(N)-m_{c}(\infty)$ vs $N$, where $m_{c}(N)\equiv\langle
m_{c}^{(i)}(N)\rangle_{i}$. By using the relation $m_{c}(N)-m_{c}(\infty)\sim
N^{-\beta_{m}/\bar{\nu}_{m}}$, $\beta_{m}/\bar{\nu}_{m}\approx 0.05$ for
$h=0.2$. The error bars represent the standard deviation of $m_{c}(N)$.

The value of $t_{c}(N)$ is determined as the mean $\langle
t_{c}^{(i)}(N)\rangle_{i}$ over the different samples; it has a relatively
large error bar. The standard deviation of $\\{t_{c}^{(i)}(N)\\}$ is defined
as:

$\displaystyle\sigma_{c}^{2}(t_{c}(N))\equiv\frac{1}{N_{\rm
config}}\sum_{i=1}^{N_{\rm config}}(t_{c}^{(i)}(N)-t_{c}(N))^{2}.\quad$ (7)

We find that as $N\to\infty$, $\sigma_{c}(t_{c}(\infty))$ does not reduce to
zero, but converges to a finite value. The difference scales as
$\sigma_{c}(t_{c}(N))-\sigma_{c}(t_{c}(\infty))\sim N^{-1.02}$, as shown in
Fig. 7(a). This corresponds to $\sigma_{c}(m_{c}(N))-\sigma(m_{c}(\infty))\sim
N^{-0.09}$ with a non-zero $\sigma_{c}(m_{c}(\infty))$. However, for the
$r$-ER model, $\sigma_{c}(t_{c}(\infty))=0$ and $\sigma_{c}(t_{c}(N))\sim
N^{-0.5}$, as shown in Fig. 7(b).

Figure 7:  (a) Plot of the finite size scaling behavior of
$\sigma_{c}(t_{c}(N))-\sigma_{c}(t_{c}(\infty))$ vs $N$ for the $m$-BFW model
with $h=0.2$. Here $\sigma_{c}(t_{c}(\infty))\approx 9.4\times 10^{-4}$ is
chosen to obtain good scaling. The finite size scaling exponent is found to be
$1.02$. (b) Plot of the standard deviation
$\sigma_{c}(m_{c}(N))-\sigma_{c}(m_{c}(\infty))$ of $m_{c}(N)$ with
$\sigma_{c}(m_{c}(\infty))\approx 0.0036$. The finite-size scaling exponent is
found to be 0.09. (c) Similar plot for the $r$-ER model with $g=0.5$:
$\sigma_{c}(t_{c}(N))$ decays in a power-law manner as $\sim N^{-0.5}$.
Therefore, $\sigma_{c}(t_{c}(\infty))=0$.

Owing to the large error bars of $t_{c}(N)$ and $\sigma_{c}(t_{c}(\infty))\neq
0$, the conventional finite-size scaling theory in the region
$t>t_{c}(\infty)$ fails. Hence, we construct a novel framework for finite-size
scaling analysis in the supercritical region as follows. We plot
$(m-m_{c}(\infty))N^{\beta_{m}/\bar{\nu}_{m}}$ vs
$(t-t_{c}^{(i)}(N))N^{1/\bar{\nu}_{m}}$ for different samples and system sizes
in double logarithm scales in Fig. 8 (a). By choosing $\beta_{m}=0.059$ and
$\bar{\nu}_{m}=1.124$, we find that the data collapse reasonably well onto a
single curve. Note that the plot is based on

${\bar{t}}^{(i)}(N)\equiv t-t_{c}^{(i)}(N)$ (8)

meaning that the transition points of each sample for a given $N$ are adjusted
to a single point instead of using their mean value $t_{c}(N)$. Note that in
the conventional finite-size scaling theory, the scaling plot is drawn as
$(m-m_{c}(\infty))N^{\beta_{m}/\bar{\nu}_{m}}$ vs
$(t-t_{c}(\infty))N^{1/\bar{\nu}_{m}}$. From this plot, we find that the data
from different system sizes are not collapsed onto a single curve, as shown in
Fig. 8(b). Hence, the conventional finite-size scaling approach established
for the second-order percolation transition does not work well for the $m$-BFW
model.

Figure 8:  (a) Finite-size scaling of order parameter $m$ for $h=0.2$ using
${\bar{t}}^{(i)}(N)$. By taking $\beta_{m}=0.059$ and $\bar{\nu}_{m}=1.124$,
the data were reasonably well collapsed onto a single curve. (b) Finite-size
scaling of order parameter $m$ in the conventional method using
$t-t_{c}(\infty)$. The data did not collapse onto a single curve.

We investigate the mean cluster size of finite clusters for systems of size
$N$. For each configuration $i$, we obtain the cluster size distribution
$n_{s}^{(i)}(t)$ and then obtain the mean cluster size of the configuration
$i$,
$\bar{s}^{(i)}(t)\equiv\sum^{\prime}_{s}s^{2}n_{s}^{(i)}(t)/\sum^{\prime}_{s}sn_{s}^{(i)}(t)$.
Because the transition point $t_{c}^{(i)}(N)$ of the configuration $i$
fluctuates significantly, we first plot the mean cluster size in scaling form
as $\bar{s}^{(i)}(t)N^{-\gamma_{s}^{\prime}/\bar{\nu}_{s}^{\prime}}$ vs
$(t_{c}^{(i)}(N)-t)N^{1/\bar{\nu}_{s}^{\prime}}$. Next, we take average of
configuration index $i$, and obtain the plot shown in Fig. 9 (a) and (b),
where $\langle s\rangle\equiv\sum_{i}{\bar{s}}^{(i)}/N_{s}$. $N_{s}$
represents the number of configurations. In Fig. 9 (a) and (b), the scaling
plot $\langle s(t)\rangle N^{-\gamma_{s}/\bar{\nu}_{s}}$ vs
$|{\bar{t}}^{(i)}(N)|N^{1/\bar{\nu}_{s}}$ exhibits data collapse behavior for
different samples and system sizes, with the choice of
$\gamma_{s}^{\prime}=1.015$ and $\bar{\nu}_{s}^{\prime}=1.020$ for $t<t_{c}$,
and $\gamma_{s}=0.918$ and $\bar{\nu}_{s}=1.031$ for $t>t_{c}$ for $h=0.2$.

Figure 9:  (a) Scaling plot of mean cluster size $\langle s\rangle
N^{-\gamma_{s}^{\prime}/\bar{\nu}_{s}^{\prime}}$ vs
$(t_{c}^{(i)}(N)-t)N^{1/\bar{\nu}_{s}^{\prime}}$ for different system sizes
(a) for $t<t_{c}(\infty)$ (in the subcritical regime) using
$\gamma_{s}^{\prime}=1.015$, $\bar{\nu}_{s}^{\prime}=1.020$ and (b) for
$t>t_{c}(\infty)$ (in the supercritical regime) using $\gamma_{s}=0.918$ and
$\bar{\nu}_{s}=1.031$.  Figure 10:  (a) Scaling plot of susceptibility
$\chi_{m}N^{-\gamma_{m}^{\prime}/\bar{\nu}_{m}^{\prime}}$ vs
$(t_{c}^{(i)}(N)-t)N^{1/\bar{\nu}_{m}^{\prime}}$ for different system sizes
(a) for $t<t_{c}(\infty)$ (in the subcritical regime) using
$\gamma_{m}^{\prime}=2.041$, $\bar{\nu}_{m}^{\prime}=1.020$ and (b) for
$t>t_{c}(\infty)$ (in the supercritical regime) using $\gamma_{m}=1.034$ and
$\bar{\nu}_{m}=1.124$

Next, we examine the susceptibility representing the fluctuations of the
largest cluster $N\chi_{m}$ for $t<t_{c}$ and of the order parameter for
$t>t_{c}$, formulated as $\chi_{m}=N(\langle m^{2}\rangle-\langle
m\rangle^{2})$. In Fig. 10(a) and (b), we attempt to show the data-collapse
behavior in the scaling plot of $\chi_{m}(t)N^{-\gamma_{m}/\bar{\nu}_{m}}$ vs
$|{\bar{t}}^{(i)}(N)|N^{1/\bar{\nu}_{m}}$ with the choices of
$\gamma_{m}^{\prime}=1.020$ and $\bar{\nu}_{m}^{\prime}=1.020$ for $t<t_{c}$,
and $\gamma_{m}=1.034$ and $\bar{\nu}_{m}\approx 1.124$ for $t>t_{c}$ and
$h=0.2$. However, the data-collapse is not satisfactory, particularly somewhat
in broad range of $t>t_{c}$ near $t_{c}^{+}$. Therefore, the critical exponent
$\gamma_{m}$ for $t>t_{c}$ is not well defined. The numerical values of the
critical exponents for other $h$ parameters are listed in Table. 2.

We have verified the scaling relations between the two sets of critical
exponents for $t>t_{c}$: one set is
$\\{\tau_{s},\sigma_{s},\beta_{s},\gamma_{s},\bar{\nu}_{s}\\}$, and the other
set is $\\{\beta_{m},\gamma_{m},\bar{\nu}_{m}\\}$. We first check the scaling
relations among the sets associated with the cluster size distribution. Using
the directly measured values $\tau_{s}=2.04\pm 0.02$ in Fig. 5 and
$\gamma_{s}=0.87\pm 0.07$ in Fig. 9, we estimate
$\sigma_{s}=(3-\tau_{s})/\gamma_{s}=1.11\pm 0.11$,
$\beta_{s}=(\tau_{s}-2)/\sigma_{s}=0.04\pm 0.02$, and
$\bar{\nu}_{s}=2\beta_{s}+\gamma_{s}=0.96\pm 0.08$. We obtained
$\bar{\nu}_{s}\approx 1.031$ as a direct estimate using the data-collapse
method in Fig. 9. Thus, the two values of $\bar{\nu}_{s}$ are in agreement
within the error range.

Next, we confirm that for $t>t_{c}$, $\sigma_{s}\approx 1$. This is consistent
with the theoretical prediction $\sigma_{s}=1$, resulting from the argument
that the size distribution of finite clusters is analytic da Costa  _et al._
(2015); Cho  _et al._ (2016)(see Appendix E).

$h$ |  $\tau_{s}$ |  $\sigma_{s}$ |  $\beta_{s}$ |  $\bar{\nu}_{s}$ |  $\gamma_{s}$ |  $\gamma_{s}^{\prime}$
---|---|---|---|---|---|---
0.2 |  $2.04\pm 0.02$ |  $1.11\pm 0.11$ |  $0.04\pm 0.02$ |  $0.96\pm 0.08$ |  $0.87\pm 0.07$ |  $1.02\pm 0.02$
0.5 |  $2.21\pm 0.03$ |  $1.04\pm 0.08$ |  $0.21\pm 0.05$ |  $1.26\pm 0.19$ |  $0.76\pm 0.03$ |  $1.08\pm 0.06$
0.8 |  $2.30\pm 0.05$ |  $1.02\pm 0.11$ |  $0.30\pm 0.08$ |  $1.33\pm 0.24$ |  $0.69\pm 0.03$ |  $1.15\pm 0.03$
$h$ |  $\beta_{m}$ |  $\bar{\nu}_{m}$ |  $\gamma_{m}$ |  $\gamma_{m}^{\prime}$ |  $\zeta$ |  $\zeta^{\prime}$
0.2 |  $0.06\pm 0.01$ |  $1.18\pm 0.06$ |  $1.05\pm 0.07$ |  $2.08\pm 0.07$ |  $1.01\pm 0.01$ |  $1.99\pm 0.08$
0.5 |  $0.24\pm 0.01$ |  $1.42\pm 0.11$ |  $0.93\pm 0.12$ |  $2.16\pm 0.02$ |  $1.02\pm 0.04$ |  $1.79\pm 0.03$
0.8 |  $0.33\pm 0.02$ |  $1.51\pm 0.13$ |  $0.85\pm 0.16$ |  $2.36\pm 0.05$ |  $1.09\pm 0.06$ |  $1.74\pm 0.08$
Table 2:  List of estimated exponent values for different values of $h$.
$t_{c}$ is the transition point in the thermodynamic limit estimated using
finite-size scaling (FSS) analysis. $\tau_{s}$ is the exponent of the cluster
size distribution at $t_{c}$. $\beta_{m}$ is the critical exponent of the
order parameter estimated using FSS analysis. $\zeta$ and $\zeta^{\prime}$ are
the exponents of the inter-cluster coalescence time distribution $P_{a}(z)$
estimated in the regimes $[0,t_{a}]$ and $[t_{a},t_{c}]$, respectively.

A scaling relation among the exponents
$\\{\beta_{m},\gamma_{m},\bar{\nu}_{m}\\}$;
$\bar{\nu}_{m}=2\beta_{m}+\gamma_{m}$ can also be verified. We estimate
$\beta_{m}=0.06\pm 0.01$ in Fig. 8(a), $\gamma_{m}=1.00\pm 0.10$ in Fig. 10.
Thus, we obtain $\bar{\nu}_{m}=1.12\pm 0.12$ from the above scaling relation.
This value is in agreement with the directly measured value
$\bar{\nu}_{m}\approx 1.124$ using the data-collapse method in both Fig. 8(a)
and Fig. 10(b).

Finally, we check whether the scaling relation Eq. (3), i.e.,
$\gamma_{s}+\beta_{m}=1$ holds. The origin of the scaling relation is similar
to the case of HPT with cluster pruning processes: The avalanche size there is
replaced in cluster merging processes by the mean cluster size. This relation
bridges the two sets of critical exponents. Using the values
$\gamma_{s}=0.87\pm 0.07$ and $\beta_{m}=0.06\pm 0.01$, we confirm that the
scaling relation $\gamma_{s}+\beta_{m}=1$ holds within the numerical accuracy.

As stated in subsection III.1 for the $r$-ER model, the inter-event time
(denoted as $z_{i}$) is defined as the time taken for a node to move from one
set to the other set. The inter-event time distribution $P_{I}(z)$ exhibits a
power-law decay with exponent $\zeta$ in the interval $[0,t_{a}]$ and exponent
$\zeta^{\prime}$ in the interval $[t_{a},t_{g}]$. This inter-event time can be
defined because the $r$-ER system is composed of two partitions. However, for
the $m$-BFW model, in which the system does not contain any partitions, the
inter-event time cannot be defined. Here, we use the age distribution as an
alternative. The age (also denoted as $z$) is defined as the duration of the
time during which the cluster of a given node is not linked with any other
cluster. For the $r$-ER model, the age distribution $P_{a}(z)$ behaves
similarly to $P_{I}(z)$: $P_{a}(z)\sim z^{-\zeta}$ for $t\in[0,t_{a}]$ and
$\sim z^{-\zeta^{\prime}}$ for $t\in[t_{a},t_{g}]$. We have found for the
$r$-ER model that Park  _et al._ (2019) $\zeta\approx 1.04$ and
$\zeta^{\prime}\approx 1.87$ for $g=0.2$. Recall that in $t\in[t_{a},t_{g}]$,
$s^{*}$ is constant for $t\to t_{c}(\infty)$ and the scaling exponent of the
characteristic size of finite clusters $\sigma_{s}^{\prime}=0$. Hence, the
relation $\zeta^{\prime}=4-(\tau_{s}+\sigma_{s}^{\prime})$ is reduced into
$\zeta^{\prime}=4-\tau_{s}$ in $[t_{a},t_{g}]$.

Figure 11:  Plot of the age distributions $P_{a}(z)$ in the intervals (a)
$t\in[0,t_{a}]$ and (b) $t\in[t_{a},t_{b}]$. The data are obtained from
$10^{5}$ configurations of system size $1.024\times 10^{7}$ for $h=0.2$

This result is also valid for the $m$-BFW model. Similar to the $r$-ER model,
characteristic size of the finite clusters is constant at a very short time
interval $[t_{a},t_{b}]$. Therefore, $\sigma_{s}^{\prime}=0$ at
$[t_{a},t_{b}]$ and the relation $\zeta^{\prime}=4-\tau_{s}$ also holds. Based
on the value, $\tau_{s}=2.04\pm 0.02$ we predict that $\zeta^{\prime}=1.96\pm
0.02$ for $h=0.2$, which is in good agreement with the value that was directly
measured, $1.99\pm 0.08$. Furthermore, using $\zeta\approx 1.01\pm 0.01$ shown
in Fig. 11 and $\tau_{s}\approx 2.04\pm 0.02$, we obtain
$\sigma_{s}^{\prime}=0.95\pm 0.03$ from the scaling relation
$\zeta=4-(\tau_{s}+\sigma_{s}^{\prime})$ in the interval $[0,t_{a}]$.

##  IV Conclusions

In this paper, we considered two models, called the modified $r$-ER model and
the $m$-BFW model, that exhibit hybrid percolation transitions (HPTs) in
cluster merging dynamics. Even though the dynamic rules of the cluster merging
process are different for the two models, we showed that their underlying
mechanisms are universal. An HPT with cluster merging occurs following three
regimes: In the early regime of $t\in[0,t_{a}]$, medium-sized clusters
accumulate and create a bump in cluster size distribution due to the
suppression of the growth of a giant cluster. This is consistent with a
“powder keg” Friedman and Landsberg (2009). Here, $t_{a}$ is defined as the
time at which the bump size reaches its maximum. In the intermediate regime
$t\in[t_{a},t_{g}]$ for the modified $r$-ER model and $t\in[t_{a},t_{b}]$ for
the $m$-BFW model, the medium-sized clusters merge and a giant cluster rapidly
increases. Time $t_{g}$ or $t_{b}$ is a crossover point at which the
suppression effect disappears. This crossover point is determined by the
dynamic rule itself in the modified $r$-ER model and in a self-organized
manner in $m$-BFW model. In the final regime $t\in[t_{g},t_{c}]$ or
$t\in[t_{b},t_{c}]$, the cluster-merging kinetics becomes effectively the ER
dynamics but with the initial configuration of cluster size distribution at
$t_{g}$ or $t_{b}$ determined by the evolutionary dynamics. At $t_{c}$, the
size distribution of finite clusters follows a power law and the system
exhibits critical behavior.

The critical behavior of the HPT is characterized by two sets of critical
exponents. One set comprising $\\{\beta_{m},\gamma_{m},\bar{\nu}_{m}\\}$ is
associated with a giant cluster, and the other set comprising
$\\{\tau_{s},\sigma_{s},\beta_{s},\gamma_{s},\bar{\nu}_{s}\\}$ is associated
with the size distribution of finite clusters. Whereas the two sets of
critical exponents are reduced to one set in the second-order, e.g., Bernoulli
percolation transition, they remain different for HPT, for example,
$\gamma_{m}\neq\gamma_{s}$. Unfortunately, due to the relatively large errors,
the two sets of exponents are numerically not that different. As the model
parameter $g$ and $h$ of the modified $r$-ER and $m$-BFW models, respectively,
change, the values of critical exponents continuously vary. However, the
exponents in each set follow scaling relations respectively and these two sets
are not independent, connected by a scaling relation $\gamma_{s}=1-\beta_{m}$.
The relation is the counterpart of that for HPT with the cluster pruning
process, with the critical exponents associated with avalanche size
distribution (cluster pruning) replaced by the exponents associated with the
cluster size distribution (cluster merging). This way a unified picture for
HPT is established.

The critical exponents are also related to the distribution of age, defined as
the duration of time during which the cluster of a given node is not linked
with any other cluster. The distributions accumulated in two time intervals,
$[0,t_{a}]$ and $[t_{a},t_{g}(t_{b})]$, exhibit power-law decays with
exponents $\zeta$ and $\zeta^{\prime}$, respectively. The two exponents can be
determined in terms of the exponents $\\{\tau_{s},\sigma_{s}$} and
$\sigma^{\prime}(=0)$, which is also true for the values of the remaining
critical exponents.

It is noteworthy that the transition point of $m$-BFW model is subject to
intrinsic uncertainty thus the conventional finite-size scaling theory in the
super-critical region breaks down to some extent. Therefore, we constructed a
novel framework for finite-size scaling analysis which adjusts the transition
points of each sample for each system size and identify the value
characteristic for the system as an average which obeys the scaling laws.

Systems exhibiting HPT constitute a very important class of dynamic
percolation problems from epidemic problems to interdependent networks, where
the jump in the order parameter represents a dramatic, explosive change of the
system. The present study completes the description of the HPT as we have now
a theoretical framework for both cluster pruning and cluster merging processes
leading to such a transition. In both cases we have two sets of exponents,
which are coupled by the scaling relation $\gamma_{a}+\beta_{m}=1$ and its
counterpart $\gamma_{s}+\beta_{m}=1$. We think that the universal mechanism
and features we uncovered in this paper can be used for understanding other
discontinuous transitions in complex systems such as in synchronization and
jamming transitions.

###### Acknowledgements.

This work was supported by the National Research Foundation of Korea by Grant
Nos. NRF-2020R1F1A1061326 (YSC), NRF-2014R1A3A2069005 (BK), KENTECH Research
Grant No. KRG2021-01-007, and US Army Research Office Grant No.
W911NF-23-1-0087 (RMD), EU Horizon 2020 grant agreement ERC No
810115-DYNASNET.

##  Appendix A Critical exponents and scaling relations

We present two sets of critical exponents to characterize the critical
behaviors of HPT on random complex networks. One set comprising
$\\{\beta_{m},\gamma_{m},\bar{\nu}_{m}\\}$ is associated with the giant
cluster size per node denoted with index $m$. The other set comprising
$\\{\tau_{s},\sigma_{s},\beta_{s},\gamma_{s},\bar{\nu}_{s}\\}$ is associated
with the size distribution of finite clusters denoted with index $s$. The
prime on the exponents represents the exponents in the subcritical regime
($t<t_{c}$) except for $\zeta^{\prime}$, specified in the Table below. The
definitions of critical exponents are listed in Table. A1.

Exponent |  Quantity |  Definition |  Regime
---|---|---|---
$\beta_{m}$ |  Giant cluster size per node |  $m-m_{c}\sim(t-t_{c})^{\beta_{m}}$ |  $t>t_{c}$
$\gamma_{m}^{\prime}$ |  Susceptibility |  $\chi_{m}\sim(t_{c}-t)^{-\gamma_{m}^{\prime}}$ |  $t<t_{c}$
$\gamma_{m}$ |  Susceptibility |  $\chi_{m}\sim(t-t_{c})^{-\gamma_{m}}$ |  $t>t_{c}$
$\bar{\nu}_{m}^{\prime}$ |  Correlation size |  $\xi_{s}\equiv\xi_{\ell}^{d}\sim(t_{c}-t)^{-d\nu_{m}^{\prime}}\sim(t_{c}-t)^{-\bar{\nu}_{m}^{\prime}}$ |  $t<t_{c}$
$\bar{\nu}_{m}$ |  Correlation size |  $\xi_{s}\equiv\xi_{\ell}^{d}\sim(t-t_{c})^{-d\nu_{m}}\sim(t-t_{c})^{-\bar{\nu}_{m}}$ |  $t>t_{c}$
$\tau_{s}$ |  Cluster size distribution at $t_{c}$ |  $n_{s}(t_{c})\sim s^{-\tau_{s}}$ |  At $t_{c}$
$\sigma_{s}^{\prime}$ |  Characteristic size |  $s^{*}\sim(t_{c}-t)^{-1/\sigma_{s}^{\prime}}$ |  $t<t_{c}$
$\sigma_{s}$ |  characteristic size |  $s^{*}\sim(t-t_{c})^{-1/\sigma_{s}}$ |  $t>t_{c}$
$\beta_{s}$ |  Giant cluster size per node |  $\beta_{s}=(\tau_{s}-2)/\sigma_{s}$ |  $t>t_{c}$
$\gamma_{s}^{\prime}$ |  Mean cluster size |  $\langle s\rangle\sim(t_{c}-t)^{-1/\gamma_{s}^{\prime}}$ |  $t<t_{c}$
$\gamma_{s}$ |  Mean cluster size |  $\langle s\rangle\sim(t-t_{c})^{-1/\gamma_{s}}$ |  $t>t_{c}$
$\bar{\nu}_{s}^{\prime}$ |  Correlation size |  $s^{*}\sim N^{1/\sigma_{s}\bar{\nu}_{s}^{\prime}}$ |  $t<t_{c}$
$\bar{\nu}_{s}$ |  Correlation size |  $s^{*}\sim N^{1/\sigma_{s}\bar{\nu}_{s}}$ |  $t>t_{c}$
$\zeta$ |  Inter-coalescence time distribution |  $P_{a}(z)\sim z^{-\zeta}$ |  $t<t_{a}$
$\zeta^{\prime}$ |  Inter-coalescence time distribution |  $P_{a}(z)\sim z^{-\zeta^{\prime}}$ |  $t_{a}<t<t_{g}~{}{\rm or}~{}t_{b}$
Table A1:  $\chi_{m}\equiv N(\langle m^{2}\rangle-\langle m\rangle^{2})$. The
distribution of finite clusters near $t_{c}$ is expressed as $n_{s}(t)\sim
s^{-\tau_{s}}f(s/s^{*})$ near $t_{c}$, where $f(x)$ is a scaling function with
$f(x)$ being constant for $x\ll 1$ and $f(x)\sim x^{\tau_{s}}$ for $x\gg 1$.
The characteristic size of finite clusters
$s^{*}\sim(t_{c}-t)^{-1/\sigma_{s}^{\prime}}$. Note that $\sigma_{s}^{\prime}$
has different values in $[0,t_{a}]$ and $[t_{a},t_{g}]$ or $[t_{a},t_{b}]$.
$\beta_{s}$ cannot be directly measured. The mean cluster size of finite
clusters $\langle
s\rangle\equiv\sum_{s}^{\prime}s^{2}n_{s}/\sum_{s}^{\prime}sn_{s}$ where the
prime in the summation denotes the sum over finite clusters.

The values of critical exponents vary depending on the model parameter $g$
($r$-ER model) or $h$ ($m$-BFW model), representing the suppression strength
in the way $1-g$ or $1-h$. The scaling relations among the critical exponents
hold independently of the model parameters. The scaling relations are listed
in Table. A2.

Giant cluster |  $\bar{\nu}_{m}=2\beta_{m}+\gamma_{m}$ |  |
---|---|---|---
Finite clusters |  $\bar{\nu}_{s}=2\beta_{s}+\gamma_{s}$ |  $\beta_{s}=(\tau_{s}-2)/\sigma_{s}$ |  $\gamma_{s}=(3-\tau_{s})/\sigma_{s}$
Scaling relation between two exponent sets |  $\gamma_{s}+\beta_{m}=1$ |  |
Scaling relation between the exponents of inter-coalescence time and cluster size distributions |  $\zeta^{\prime}=4-(\tau_{s}+\sigma_{s}^{\prime})$ |  $\zeta=4-\tau_{s}$ |
Table A2:  Scaling relations among the critical exponents.

##  Appendix B Difference between the original $r$-ER and the modified $r$-ER
models

Here, we explain the difference between the original and modified $r$-ER
models in their dynamic rules. Then, why the modification is necessary to make
an HPT occur is discussed.

The original $r$-ER model was proposed to study a discontinuous percolation
transition in Ref. Panagiotou  _et al._ (2011). In the model, $g$ fraction of
nodes that are contained in the smallest clusters are assigned to set $A$ and
the remaining fraction of nodes are assigned to set $B$. As clusters are
merged, clusters can be reassigned their sets. In the original model, the
number of nodes assigned to set $A$ is fixed as $\lfloor gN\rfloor$. As a
result, nodes of the cluster on the border are assigned to either set $A$ or
$B$. This assignment rule differs from that of the modified $r$-ER model,
where nodes in the cluster on the border of the two sets $A$ and $B$ are
regarded as the elements of set $A$. As a result, when $t$ reaches $t_{g}$,
then $m(t_{g})=1-g$, that is, the fractional size of the largest cluster at
time $t_{g}$ becomes the capacity $1-g$ of set $B$, the largest cluster is
regarded as the element of set $A$. Then, the partitions become unified, and
the dynamics of cluster merging is governed by the ER dynamics when $t\geq
t_{c}$.

Figure A1:  Schematic plots of the original $r$-ER model (a) and the modified
$r$-ER model (b). Circles represent clusters. Circles contained in set $A$
($B$) are marked in red (blue).

##  Appendix C Evolution of the order parameter and the cluster size
distribution in the three regimes of the modified $r$-ER model

Figure A2:  Plots of the cluster size distributions of the modified $r$-ER
model at (a) $t\ll t_{a}$, (b) $t=t_{a}$, (c) $t=t_{c}$, and (d) $t>t_{c}$.
Inset: Schematic plots of the order parameter $m(t)$ vs $t\equiv L/N$ with a
point ($\star$) indicating the corresponding time. (e) Plot of the order
parameter with characteristic points $\\{t_{a},t_{g},t_{c}\\}$ in the
evolution of a single sample of the modified $r$-ER model with a system size
$N=1.024\times 10^{4}$ for $g=0.5$. Similar plots were presented in Refs.Cho
_et al._ (2016); Park  _et al._ (2019)

We plot the evolution of the order parameter and the cluster size distribution
of the modified $r$-ER model throughout the three regimes in Fig. A2. This
plot is shown for the comparison to Fig. 2 of the $m$-BFW model.

##  Appendix D Difference between the original and $m$-BFW models

Figure A3:  (a) Plot of the acceptance rate $L/u$ (left axis) and the limit
$k$ of the largest cluster size (right axis) vs time $t$ for the original BFW
model. $h(k)$ decreases regularly. Thus, the limit $k$ increases regularly,
and the step widths and heights respectively decrease and increase regularly,
respectively. The acceptance rate, i.e., $L/u$, also decreases. (b) Plots of
the age distribution as a function of age $z$. The age distribution for
$t\in[t_{a},t_{b}]$ does not exhibit power-law behavior, and thus the
criticality is not formed for $t>t_{c}$ and the percolation transition is
discontinuous.

The main difference between the original and $m$-BFW models lies in the lower
bound $h$ of the acceptance rate $L/u$. For the $m$-BFW model, $h$ is
constant, whereas for the original BFW model, $h(k)$ depends on $k$ as
$h(k)=1/2+(2k)^{-1/2}$. Accordingly, the acceptance rate $L/u$ behaves
differently as a function of time $t$ as shown in Fig. A3. We found that
$h(k)$ decreases regularly, and thus the upper bound of cluster size denoted
as $k$ increases regularly; the step widths and heights respectively decrease
and increase regularly. The acceptance rate, i.e., $L/u$, also decreases.
Accordingly, the age distribution does not show a power-law behavior in the
interval $t\in[t_{a},t_{b}]$.

##  Appendix E Derivation of $\sigma_{s}=1$

Let $p_{s}(\hat{t})$ be the probability that a randomly selected node is an
element of a cluster of size $s$ at relative time $\hat{t}\equiv t-t_{c}$.
Then $p_{s}(0)\sim s^{1-\tau}$ and $\sum_{s=1}^{\prime}p_{s}(0)=1-m_{c}$ for
HPT. For $\hat{t}\geq 0$, the kinetics of the cluster merging of the modified
$r$-ER model is governed by the ER dynamics with the initial configuration
$p_{s}(0)$. For $\hat{t}>0$, $p_{s}(\hat{t})$ is analytic with respect to
$\hat{t}$, and thus $p_{s}(\hat{t})$ can be written as

$p_{s}(\hat{t})=B_{0}(s)+B_{1}(s)\hat{t}+B_{2}(s)\hat{t}^{2}+\cdots,$ (9)

where $B_{0}(s)=p_{s}(0)$. Using the property
$\sum_{s=1}^{\prime}B_{0}(0)=1-m_{c}$, one can find that $B_{1}(s)\sim
s^{2-\tau}$ and $B_{1}(s)\sim s^{3-\tau}$, etc. It is straightforward to
derive $B_{n}(s)\sim s^{n+1-\tau}$ for general $n\geq 0$ using the relation
$B_{n+1}(s)\sim sB_{n}(s)$ when $s$ is large. Therefore, $p_{s}(\hat{t})$ can
be written as $p_{s}(\hat{t})=s^{1-\tau}\sum_{n\geq 0}b_{n}(s\hat{t})^{n}$
with constant $b_{n}$. The derivation of this result was presented in detail
in the supplementary information of Ref. Cho  _et al._ (2016). Therefore,
$p_{s}(\hat{t})$ can be written in a scaling form,
$p_{s}(\hat{t})=s^{1-\tau}f(s\hat{t}^{1/\sigma_{s}})$ with $\sigma_{s}=1$.
This approach can also be used in the $m$-BFW model.

## References

  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Stanley (1971) H.E. Stanley,  _Introduction to Phase Transitions and Critical Phenomena_, The International Series of Monographs on Physics (Oxford University Press, 1971).
  * Blossey and Indekeu (1995) R. Blossey and J. O. Indekeu, “Diverging length scales at first-order wetting transitions,” Phys. Rev. E 52, 1223–1226 (1995).
  * Kafri  _et al._ (2000) Yariv Kafri, David Mukamel, and Luca Peliti, “Why is the dna denaturation transition first order?” Phys. Rev. Lett. 85, 4988–4991 (2000).
  * Bar and Mukamel (2014) Amir Bar and David Mukamel, “Mixed-order phase transition in a one-dimensional model,” Phys. Rev. Lett. 112, 015701 (2014).
  * Thouless (1969) D. J. Thouless, “Long-range order in one-dimensional ising systems,” Phys. Rev. 187, 732–733 (1969).
  * Dyson (1971) Freeman J Dyson, “An ising ferromagnet with discontinuous long-range order,” Commun. Math. Phys. 21, 269–283 (1971).
  * Aizenman  _et al._ (1988) Michael Aizenman, JT Chayes, L Chayes, and CM Newman, “Discontinuity of the magnetization in one-dimensional 1/¦x- y¦ 2 ising and potts models,” J. Stat. Phys. 50, 1–40 (1988).
  * Jang  _et al._ (2015) S. Jang, J. S. Lee, S. Hwang, and B. Kahng, “Ashkin-teller model and diverse opinion phase transitions on multiplex networks,” Phys. Rev. E 92, 022110 (2015).
  * Toninelli  _et al._ (2006) Cristina Toninelli, Giulio Biroli, and Daniel S. Fisher, “Jamming percolation and glass transitions in lattice models,” Phys. Rev. Lett. 96, 035702 (2006).
  * Janssen and Stenull (2016) Hans-Karl Janssen and Olaf Stenull, “First-order phase transitions in outbreaks of co-infectious diseases and the extended general epidemic process,” EPL 113, 26005 (2016).
  * Choi  _et al._ (2017a) Wonjun Choi, Deokjae Lee, and B. Kahng, “Mixed-order phase transition in a two-step contagion model with a single infectious seed,” Phys. Rev. E 95, 022304 (2017a).
  * Hasegawa and Nemoto (2018) Takehisa Hasegawa and Koji Nemoto, “Sudden spreading of infections in an epidemic model with a finite seed fraction,” Eur. Phys. J. B 91, 1–8 (2018).
  * Pazó (2005) Diego Pazó, “Thermodynamic limit of the first-order phase transition in the kuramoto model,” Phys. Rev. E 72, 046211 (2005).
  * Gómez-Gardeñes  _et al._ (2011) Jesús Gómez-Gardeñes, Sergio Gómez, Alex Arenas, and Yamir Moreno, “Explosive synchronization transitions in scale-free networks,” Phys. Rev. Lett. 106, 128701 (2011).
  * Coutinho  _et al._ (2013) B. C. Coutinho, A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes, “Kuramoto model with frequency-degree correlations on complex networks,” Phys. Rev. E 87, 032106 (2013).
  * Park and Kahng (2019) Jinha Park and B. Kahng, “Abnormal hybrid phase transition in the passively competing kuramoto model,” Phys. D 399, 186–192 (2019).
  * Basnarkov and Urumov (2008) Lasko Basnarkov and Viktor Urumov, “Kuramoto model with asymmetric distribution of natural frequencies,” Phys. Rev. E 78, 011113 (2008).
  * Park  _et al._ (2020) Jinha Park, Sudo Yi, and B. Kahng, “Hysteresis and criticality in hybrid percolation transitions,” Chaos 30, 051102 (2020).
  * Panagiotou  _et al._ (2011) Konstantinos Panagiotou, Reto Spöhel, Angelika Steger, and Henning Thomas, “Explosive percolation in erdős-rényi-like random graph processes,” Electron. Notes Discrete Math. 38, 699–704 (2011).
  * Buldyrev  _et al._ (2010) Sergey V Buldyrev, Roni Parshani, Gerald Paul, H Eugene Stanley, and Shlomo Havlin, “Catastrophic cascade of failures in interdependent networks,” Nature 464, 1025–1028 (2010).
  * Baxter  _et al._ (2012) G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “Avalanche collapse of interdependent networks,” Phys. Rev. Lett. 109, 248701 (2012).
  * Dorogovtsev  _et al._ (2006) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “$k$-core organization of complex networks,” Phys. Rev. Lett. 96, 040601 (2006).
  * Choi  _et al._ (2021) K. Choi, Wonjun Choi, and B. Kahng, “A hybrid percolation transition at a finite transition point in scale-free networks,” Chaos 31, 053128 (2021).
  * Zhou  _et al._ (2014a) Dong Zhou, Amir Bashan, Reuven Cohen, Yehiel Berezin, Nadav Shnerb, and Shlomo Havlin, “Simultaneous first- and second-order percolation transitions in interdependent networks,” Phys. Rev. E 90, 012803 (2014a).
  * Schwarz  _et al._ (2006) J. M. Schwarz, A. J. Liu, and L. Q. Chayes, “The onset of jamming as the sudden emergence of an infinite k-core cluster,” Europhys. Lett. 73, 560 (2006).
  * Son  _et al._ (2012) Seung-Woo Son, Golnoosh Bizhani, Claire Christensen, Peter Grassberger, and Maya Paczuski, “Percolation theory on interdependent networks based on epidemic spreading,” Europhys. Lett. 97, 16006 (2012).
  * Choi  _et al._ (2017b) K. Choi, Deokjae Lee, Y. S. Cho, J. C. Thiele, H. J. Herrmann, and B. Kahng, “Critical phenomena of a hybrid phase transition in cluster merging dynamics,” Phys. Rev. E 96, 042148 (2017b).
  * Alert  _et al._ (2017) Ricard Alert, Pietro Tierno, and Jaume Casademunt, “Mixed-order phase transition in a colloidal crystal,” PNAS 114, 12906–12909 (2017).
  * Araújo  _et al._ (2014) N Araújo, Peter Grassberger, Byungnam Kahng, KJ Schrenk, and Robert M Ziff, “Recent advances and open challenges in percolation,” Eur. Phys. J. 223, 2307–2321 (2014).
  * Dorogovtsev  _et al._ (2008) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “Critical phenomena in complex networks,” Rev. Mod. Phys. 80, 1275–1335 (2008).
  * Lee  _et al._ (2018) Deokjae Lee, B Kahng, YS Cho, K-I Goh, and D-S Lee, “Recent advances of percolation theory in complex networks,” Korean Phys. Soc. 73, 152–164 (2018).
  * Erdős and Rényi (1960) Paul Erdős and Alfréd Rényi, “On the evolution of random graphs,” Publ. Math. Inst. Hung. Acad. Sci 5, 17–60 (1960).
  * Achlioptas  _et al._ (2009) Dimitris Achlioptas, Raissa M. D’Souza, and Joel Spencer, “Explosive percolation in random networks,” Science 323, 1453–1455 (2009).
  * Bollobás (2001) Béla Bollobás,  _Random Graphs_, 2nd ed., Cambridge Studies in Advanced Mathematics (Cambridge University Press, 2001).
  * Stauffer and Aharony (1992) Dietrich Stauffer and Ammon Aharony,  _Introduction to percolation theory, 2nd ed_ (Taylor and Francis, 1992).
  * Christensen and Moloney (2005) Kim Christensen and Nicholas R Moloney,  _Complexity and Criticality_ (Imperial College Press, 2005).
  * D’Souza and Nagler (2015) Raissa M D’Souza and Jan Nagler, “Anomalous critical and supercritical phenomena in explosive percolation,” Nat. Phys 11, 531–538 (2015).
  * Ben-Naim and Krapivsky (2005) E. Ben-Naim and P. L. Krapivsky, “Kinetic theory of random graphs: From paths to cycles,” Phys. Rev. E 71, 026129 (2005).
  * Lee  _et al._ (2017) Deokjae Lee, Wonjun Choi, Janos Kertész, and Byungnam Kahng, “Universal mechanism for hybrid percolation transitions,” Sci. Rep. 7, 1–7 (2017).
  * Zhou  _et al._ (2014b) Dong Zhou, Amir Bashan, Reuven Cohen, Yehiel Berezin, Nadav Shnerb, and Shlomo Havlin, “Simultaneous first- and second-order percolation transitions in interdependent networks,” Phys. Rev. E 90, 012803 (2014b).
  * Lee  _et al._ (2016a) Deokjae Lee, S. Choi, M. Stippinger, J. Kertész, and B. Kahng, “Hybrid phase transition into an absorbing state: Percolation and avalanches,” Phys. Rev. E 93, 042109 (2016a).
  * Chalupa  _et al._ (1979) J Chalupa, P L Leath, and G R Reich, “Bootstrap percolation on a bethe lattice,” J. Phys. C 12, L31–L35 (1979).
  * Baxter  _et al._ (2015) G. J. Baxter, S. N. Dorogovtsev, K.-E. Lee, J. F. F. Mendes, and A. V. Goltsev, “Critical dynamics of the $k$-core pruning process,” Phys. Rev. X 5, 031017 (2015).
  * Lee  _et al._ (2016b) Deokjae Lee, Minjae Jo, and B. Kahng, “Critical behavior of $k$-core percolation: Numerical studies,” Phys. Rev. E 94, 062307 (2016b).
  * Cho  _et al._ (2016) Y. S. Cho, J. S. Lee, H. J. Herrmann, and B. Kahng, “Hybrid percolation transition in cluster merging processes: Continuously varying exponents,” Phys. Rev. Lett. 116, 025701 (2016).
  * Park  _et al._ (2019) Jinha Park, Sudo Yi, K. Choi, Deokjae Lee, and B. Kahng, “Interevent time distribution, burst, and hybrid percolation transition,” Chaos 29, 091102 (2019).
  * Friedman and Landsberg (2009) Eric J. Friedman and Adam S. Landsberg, “Construction and analysis of random networks with explosive percolation,” Phys. Rev. Lett. 103, 255701 (2009).
  * Cho  _et al._ (2010) Y. S. Cho, B. Kahng, and D. Kim, “Cluster aggregation model for discontinuous percolation transitions,” Phys. Rev. E 81, 030103 (2010).
  * Bohman  _et al._ (2004) Tom Bohman, Alan Frieze, and Nicholas C. Wormald, “Avoidance of a giant component in half the edge set of a random graph,” Random Struct. Algorithms 25, 432–449 (2004).
  * Chen and D’Souza (2011) Wei Chen and Raissa M. D’Souza, “Explosive percolation with multiple giant components,” Phys. Rev. Lett. 106, 115701 (2011).
  * da Costa  _et al._ (2015) R. A. da Costa, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “Solution of the explosive percolation quest. ii. infinite-order transition produced by the initial distributions of clusters,” Phys. Rev. E 91, 032140 (2015).