Self-similar dynamics in percolation and sandpile

Abstract

Spatial self-similarity is a hallmark of critical phenomena. We study the dynamic process of percolation, in which bonds are incrementally added to an initially empty lattice until the system becomes fully occupied. By tracking the gap – the size increment of clusters upon bond addition – and the corresponding merged cluster, we identify scale-invariant temporal patterns in both quantities throughout a large portion of the process. This reveals a form of temporal self-similarity that has not been reported before. We further establish quantitative relations between the dynamic scaling exponents and the conventional static critical exponents, which enable the determination of critical behavior without prior knowledge of the critical point. The same self-similar dynamics is observed in both bond and site percolation on lattices and networks, and extends to other systems such as explosive and rigidity percolation. Moreover, similar temporal scaling is found in the initial nonequilibrium evolution of the Bak−Tang−Wiesenfeld sandpile model, suggesting a dynamic critical behavior distinct from its equilibrium state. These results provide a unified framework for understanding critical dynamics and may find applications in a broad range of complex systems.

Graphical abstract

Keywords

percolation / sandpile / dynamic / self-similarity / complex systems / critical phenomena

1 Introduction

Critical phenomena describe the behaviors of systems undergoing continuous phase transitions [1], where macroscopic properties change dramatically in response to small variations of control parameters. A key feature is spatial self-similarity, where structures at different scales appear statistically identical. Percolation is a paradigmatic model for studying critical phenomena [2]. As the occupation probability p approaches the percolation threshold p_c, the correlation length diverges as \xi \sim |p-p_c{|}^{-\nu }, and self-similar spatial structures emerge up to a characteristic size s_{\xi }\sim {\xi }^{d_f}\sim |p-p_c{|}^{-1/\sigma }, with d_f the fractal dimension. The cluster number density obeys

(1)n(s)=s^{-\tau }\tilde{n}(s/s_{\xi }),

with \tau the Fisher exponent and \tilde{n}(\cdot ) a universal scaling function. Most critical behaviors can be captured by two independent exponents, such as \nu and d_f, from which others can be derived via scaling relations. In d dimensions, 1/\sigma =\nu d_f and \tau =1+d/d_f. In two dimensions (2D), \nu =4/3 and d_f=91/48 give \sigma =36/91 and \tau =187/91 [3–7]. Changing spatial dimension or cluster formation rules, such as in explosive percolation (EP) [8, 9] and rigidity percolation [10–13], can alter critical exponents, but the scaling relations hold across models.

Numerical studies of critical phenomena are generally performed on finite systems, where self-similar behavior is characterized by finite-size scaling (FSS) theory. The central assumption of FSS is that, for a system of linear size L, when the control parameter p lies within a finite-size critical window |p-p_c|\in \mathcal{O}(L^{-1/\nu }), the correlation length saturates at \mathcal{O}(L). As a consequence, singular critical observables manifest as power-law scaling with respect to the system size L. For example, the size of the largest cluster C_1, which represents the characteristic cluster size s_{\xi }, follows the scaling relation C_1\sim L^{d_f}. Such FSS relations provide a standard numerical route to determine critical exponents, provided that critical observables can be accurately measured within the finite-size critical window for a series of system sizes L.

However, this critical window \mathcal{O}(L^{-1/\nu }) shrinks rapidly as L increases. As a result, even small uncertainties in the estimated critical point p_c can in principle shift large-scale simulations outside the true scaling regime, thereby invalidating the application of FSS. This makes FSS analyses sensitive to the precise, a priori determination of p_c. Moreover, even when p_c is determined with high accuracy, strong finite-size corrections or large sample-to-sample fluctuations often obscure the asymptotic scaling behavior in many complex systems. For instance, EP was initially misidentified as a discontinuous phase transition due to its anomalous FSS behavior [8, 14–19]. Similarly, in 2D rigidity percolation, reliable estimates of critical exponents remain challenging because of pronounced finite-size effects [10–13, 20–26].

These limitations reflect a more fundamental issue: conventional FSS relies on finely tuning the control parameter to a critical point and extracting scaling behavior from essentially static snapshots. In contrast, when a system evolves dynamically and the control parameter is continuously varied, it naturally traverses the critical regime and encodes rich temporal information, which is largely discarded in static analyses. Motivated by this perspective, exploring the inherent dynamic structure of cluster formation provides a complementary route to characterize critical behavior, without requiring a priori knowledge of the critical point. Such a dynamic framework also better reflects many real-world processes, ranging from power grids [27] and climate dynamics [28] to traffic congestion [29] and information spreading [30].

In this paper, we uncover a temporal self-similarity that emerges from the dynamics of cluster merging in percolation. For bond percolation, we define evolutionary time t as the fraction of occupied bonds B/E, where B is the number of inserted bonds and E is the total number of bonds. This simple bond-insertion process, while monotonically related to the overall bond occupied probability p, reveals underlying kinetic correlations. At each step, we record the minimal size increment in clusters (termed gap) and the size of the resulting merged cluster. The temporal evolution reveals self-similar dynamics: both the gap-size distribution {\mathcal{P}}_G(s,L) and the cluster-size distribution {\mathcal{P}}_S(s,L) conform to the scaling form of Eq. (1); see Figs. 1(a) and (b). The cutoff size is s_{\xi }\sim L^{d_f}, while the Fisher exponents are {\tau }_G^{\prime }=\tau for gaps and {\tau }_S^{\prime }=\tau -1+\sigma for clusters. Dynamic observables exhibit clean FSS behaviors, e.g., the largest gap G_{max}\sim L^{d_f}.

This dynamic provides a novel, universally applicable, and highly efficient framework for extracting critical properties, crucially bypassing the need for precise prior determination of the critical point p_c. This addresses a major technical bottleneck in conventional FSS analysis. A simple example is 1D percolation, where p_c=1 precludes any critical behavior in conventional FSS analyses. However, its temporal dynamics clearly exhibits self-similarity, with both {\mathcal{P}}_G(s,L) and {\mathcal{P}}_S(s,L) following a power law characterized by {\tau }_G^{\prime }={\tau }_S^{\prime }=2, yielding \nu =d_f=1. Applied in EP, this approach shows that the system still obeys the standard FSS theory, and allows for the extraction of high-precision critical exponents. Furthermore, the dynamic analysis in rigidity percolation immediately reveals a previously unreported self-similar cascade of cluster mergings [26, 31], underscoring the transformative power of the temporal self-similarity framework.

To demonstrate the broad scope and universality of this dynamic scaling, we apply the temporal framework to the Bak−Tang−Wiesenfeld (BTW) sandpile model [32–35], a paradigmatic example of self-organized criticality (SOC). In this model, grains are added one by one to randomly selected sites of an initially empty square lattice with open boundaries. Each site can hold at most three grains, and when a fourth grain arrives, it topples, transferring one grain to each of its four neighbors. This toppling may trigger further topplings, leading to a cascade of events called an avalanche.

We focus on the early, transient dynamics leading up to the first spanning avalanche, before the system settles into the stationary SOC state. At each step, we record both the avalanche size S (the total number of topplings) and the avalanche area A (the number of distinct sites toppled). Both distributions {\mathcal{P}}_S and {\mathcal{P}}_A display clear power-law forms. The discovery that this same functional form of temporal self-similarity exists in the non-equilibrium BTW model strongly supports the hypothesis of a broad, underlying temporal universality in critical dynamics, regardless of whether the system is externally tuned or self-organized. Interestingly, while the avalanche area A is governed by a single fractal dimension, the avalanche size S does not collapse with a single exponent, suggesting a richer, multifractal behavior akin to that reported in earlier studies [36]. This temporal approach thus provides a fresh lens for understanding the complex scaling properties of SOC dynamics.

The remainder of this paper is organized as follows. In Section 2, we introduce the gap-dynamics framework and define the key observables. Sections 3 and 4 present the theoretical predictions and numerical results for various percolation systems. Section 5 applies the dynamic framework to the BTW sandpile model. Finally, Section 6 provides a brief discussion.

2 Gap dynamics and observables

2.1 Dynamic bond-addition process

We introduce a unified dynamic framework applicable to a broad class of percolation-type systems, including bond, site, explosive, and rigidity percolation on lattices and networks. The system starts from an empty configuration and evolves through the sequential addition of elementary elements, such as bonds or sites, one at a time. After each addition, the newly introduced element may connect multiple existing clusters, resulting in the formation of a larger merged cluster. The number of added elements (bonds or sites) is denoted as B. This process continues until a tunable endpoint B_e, corresponding to full occupation when B_e equals the total number of bonds E or sites V in the system.

To facilitate comparison across different systems, we define a normalized time-like variable

(2)t\equiv {\displaystyle \frac{B}{E}},\text{or}t\equiv {\displaystyle \frac{B}{V}},

depending on whether bonds or sites are added. The corresponding endpoint is denoted by t_e, with a maximum value of 1 corresponding to the fully occupied state. Numerically, the dimensionless variable t\in [0,1] is equivalent to the bond/site occupation probability p. However, they play distinct conceptual roles: while p characterizes a static configuration at a given occupation level, t parameterizes the dynamical evolution of the system through successive addition events and thus serves as a natural time-like variable in the dynamic process.

For Bernoulli bond or site percolation, an analogous dynamic process underlies the celebrated Newman−Ziff algorithm [37]. Moreover, similar dynamic constructions naturally arise in a broad class of percolation-related systems, including EP [8] and growing graphs [38, 39].

2.2 Observables

At each merging event during the dynamic process, we characterize the structural change using two quantities: the merged cluster size \mathcal{S} and the gap size \mathcal{G}. The merged cluster size

(3)\mathcal{S}={\displaystyle \sum _i{s}_i}

measures the size of the cluster formed after a merging event, where s_i denote the sizes of the clusters involved. To quantify the effective impact of the merging operation itself, we define the gap size as

(4)\mathcal{G}=\mathcal{S}-s_{\mathrm{m}\mathrm{a}\mathrm{x}},

which corresponds to the contribution from all but the largest participating cluster s_{\mathrm{m}\mathrm{a}\mathrm{x}}. By excluding the trivial contribution from a dominant cluster, \mathcal{G} isolates the genuine structural change induced by the merging event and remains sensitive to the dynamics of finite clusters even in the presence of a giant cluster. If the added element is internal, i.e., it does not merge distinct clusters, then \mathcal{G}=0 and \mathcal{S} reduces to the size of the cluster involved.

During the dynamic process, the gap size \mathcal{G}(t) typically exhibits an overall growth-decay trend: it increases at early stages, reaches a maximum value {\mathcal{G}}_{\mathrm{m}\mathrm{a}\mathrm{x}} at a characteristic dynamic pseudocritical time t_L, and then decreases, while displaying fluctuations at the level of individual merging events. The ensemble-averaged quantities G_{\mathrm{m}\mathrm{a}\mathrm{x}}\equiv ⟨{\mathcal{G}}_{\mathrm{m}\mathrm{a}\mathrm{x}}⟩, together with the mean and variance of t_L, naturally characterize the dynamic critical behavior and form the basis of the FSS analysis [40–43].

In this work, we are interested in the probability distribution of gap sizes during the dynamic process. We sample the distribution {\mathcal{P}}_G(s,L), where {\mathcal{P}}_G(s,L)\mathrm{\Delta }s denotes the normalized count of gap events with sizes in the interval (s-\mathrm{\Delta }s/2,s+\mathrm{\Delta }s/2]. To properly capture the broad range of gap sizes, the bin width \mathrm{\Delta }s is chosen to increase as a power law of s. While {\mathcal{P}}_G(s,L) can in principle be obtained from a single stochastic realization, reliable statistics generally require averaging over many independent realizations.

The average and squared-average of the gap size over an entire dynamic process are then given by the first and second moments of {\mathcal{P}}_G(s,L),

(5)G_{\mathrm{a}\mathrm{v}\mathrm{e}}(L)={\displaystyle \sum _{s}s{\mathcal{P}}_G(s,L)=⟨{\displaystyle \frac{1}{B_e}}{\displaystyle \sum _{B=1}^{B_e}\mathcal{G}(B)⟩,}}

(6)G_2(L)={\displaystyle \sum _s{s}^2{\mathcal{P}}_G(s,L)=⟨{\displaystyle \frac{1}{B_e}}{\displaystyle \sum _{B=1}^{B_e}{\mathcal{G}}^2(B)⟩,}}

where ⟨\cdot ⟩ denotes an average over independent realizations ({10}^5–{10}^6 for typical system sizes). Analogous quantities S_{\mathrm{a}\mathrm{v}\mathrm{e}}, S_2, and the dynamic distribution {\mathcal{P}}_S(s,L) are defined for the merged cluster size \mathcal{S}.

It is natural to probe how a critical configuration responds to a local bond or site addition, which provides a susceptibility-like measure of the instantaneous response of critical states to a minimal structural perturbation. At the percolation threshold t_c=p_c, we therefore consider the effect of adding a single bond or site at random to a critical configuration and measure the resulting quantities \mathcal{G} and \mathcal{S}. This single-addition operation acts as a local probe of the critical structure, and the added bond or site is not retained after the measurement. Repeating this procedure over many independent critical configurations yields a large ensemble of values of \mathcal{G} and \mathcal{S}, which characterize the instantaneous response of critical states. For notational convenience, we refer to the corresponding distributions as the static distributions and denote them by P_G(s,L) and P_S(s,L), using normal fonts to distinguish them from the dynamic distributions {\mathcal{P}}_G(s,L) and {\mathcal{P}}_S(s,L) obtained along the evolving process.

We also note that the cluster-number density n(s,L), a standard static quantity, is defined as the number of clusters of size s per unit volume V=L^d in a configuration, and should be distinguished from the static distribution P_S(s,L).

3 Self-similarity in percolation dynamics

3.1 Static critical configuration

Similar to the cluster-number density n(s,L) in Eq. (1), the static distributions at criticality should obey the FSS forms

(7)P_X(s,L)=s^{-{\tau }_X}{\tilde{P}}_X(s/L^{d_f}),

where X denotes either the gap size G or the cluster size S. Thus, the first moment of P_X(s,L) gives the scaling forms of the average sizes of cluster and gap,

(8)S_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{(2-{\tau }_S)d_f},

(9)G_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{(2-{\tau }_G)d_f}.

By definition, for a static configuration the size \mathcal{S} equivalently represents the cluster size accessed by a randomly chosen bond or site. Its ensemble average S_{\mathrm{a}\mathrm{v}\mathrm{e}} corresponds to the magnetic susceptibility of the percolation transition, i.e., the second moment of the cluster-number density n(s,L) [2], with the known FSS

(10)S_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{2d_f-d},

where d is the system dimension. Similarly, the average gap size G_{\mathrm{a}\mathrm{v}\mathrm{e}} corresponds to the typical size of the smaller cluster connected by a randomly chosen bond. Equivalently, for site percolation, it measures all clusters accessed by a site except the largest one. Its FSS form has been derived in Ref. [44],

(11)G_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{d_f-d+1/\nu },

where \nu is the correlation-length exponent. Comparing Eqs. (8) and (9) with Eqs. (10) and (11) leads to the relations between the Fisher exponents

(12){\tau }_S=\tau -1,

(13){\tau }_G=\tau -\sigma ,

where \tau =1+d/d_f is the standard Fisher exponent as in Eq. (1) and \sigma =1/(\nu d_f) is the critical exponent of the characteristic cluster size s_{\xi }. Since the probability that a randomly chosen bond or site belongs to a cluster of size s is proportional to s, the relation P_S(s,L)\sim sn(s,L) holds, recovering Eq. (12). These relations show how the static gap and cluster distributions encode the standard Fisher exponents of percolation, connecting the single-bond or -site response to the global critical scaling.

In 2D, with d_f=91/48 and \nu =4/3, one has \tau =187/91 and \sigma =36/91, yielding {\tau }_S=96/91 and {\tau }_G=151/91. Figures 1(c1) and (c2) show the simulation results for the static cluster-size and gap-size distributions P_S(s,L) and P_G(s,L) obtained from bond percolation on square lattices. The measured distributions follow the predicted power-law behaviors P_G\sim s^{-{\tau }_G} and P_S\sim s^{-{\tau }_S} over a broad range, in excellent agreement with the theoretical exponents. The insets present the FSS of their moments, confirming the expected scaling: S_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{(2-{\tau }_S)d_f}\sim L^{43/24}, G_{\mathrm{a}\mathrm{v}\mathrm{e}}\sim L^{(2-{\tau }_G)d_f}\sim L^{31/48}, S_2\sim L^{(3-{\tau }_S)d_f}\sim L^{59/16}, and G_2\sim L^{(3-{\tau }_G)d_f}\sim L^{61/24}. These results provide a clear verification of the static critical scaling and validate the Fisher exponents derived above.

3.2 Dynamic process

Figures 1(a1, a2) and (b1, b2) show the dynamic distributions {\mathcal{P}}_S(s,L) and {\mathcal{P}}_G(s,L) for bond percolation on square lattices, measured up to t_e=t_L and t_e=1, respectively. Both distributions exhibit clear algebraic scaling, indicating dynamic self-similarity. However, the associated Fisher exponents differ from their static counterparts, with {\tau }_S^{\prime }=132/91 and {\tau }_G^{\prime }=187/91. Careful analysis shows that these dynamic exponents satisfy a simple shift relation (to be derived later)

(14){\tau }_S^{\prime }={\tau }_S+\sigma =\tau +\sigma -1,

(15){\tau }_G^{\prime }={\tau }_G+\sigma =\tau .

The fact that {\tau }_G^{\prime }=\tau is highly interesting. It implies that the dynamic gap-size distribution obeys the same clean scaling as the static cluster-number density at criticality, even beyond the percolation threshold (t_e>t_L). As shown in the inset, its moments also follow the same FSS behavior as those of the static cluster-number density. This provides a powerful framework to extract critical properties without prior knowledge of the precise critical point. Such robustness originates from the temporal accumulation inherent in the dynamic process and is consistent with the general framework of critical phenomena. A more detailed explanation will be given later.

In contrast, the dynamic cluster-size distribution {\mathcal{P}}_S(s,L) exhibits a Fisher exponent {\tau }_S^{\prime }=\tau -1+\sigma =(d+1/\nu )/d_f, which encodes information from both the fractal dimension d_f and the correlation-length exponent \nu. This reflects the richer structure of temporal accumulation in {\mathcal{P}}_S, incorporating contributions from both the critical regime and the approach-to-critical stages of the dynamic process.

Furthermore, when the process is terminated in the supercritical phase (e.g., at t_e=1), {\mathcal{P}}_S(s,L) acquires an additional contribution from the giant cluster C_1\sim L^d, reflecting the nonzero probability that a bond belongs to the largest cluster in the supercritical regime. In Fig. 1(b1), this contribution manifests as a distinct bump for large s, whose height decays as \sim s^{-1}. This feature can also be interpreted within the standard framework of phase transition theory and will be derived in detail later. In contrast, the gap-size distribution {\mathcal{P}}_G(s,L) remains free from such a bump, since the definition of the gap naturally excludes any contribution from the giant cluster, see Fig. 1(b2).

3.3 Derivation of dynamic Fisher exponents

We now derive the dynamic Fisher exponent relations, Eqs. (14) and (15), within the framework of phase transition and FSS theories. The key idea is that the dynamic distributions can be regarded as the temporal accumulation of the corresponding static ones. Given an evolution time t for the dynamic process, as long as the instantaneous correlation length \xi (t) is much larger than the lattice spacing, each static observable at time t follows the critical scaling form characterized by relations such as \xi (t)\sim |t-t_c{|}^{-\nu } and P_X(s,\xi (t))\sim s^{-{\tau }_X}{\tilde{P}}_X(s/s_{\xi }(t)). Their accumulation over time thus naturally preserves the scaling behavior, leading to well-defined dynamic exponents.

3.3.1 Universal dynamic scaling of finite clusters

We begin by considering the process that terminates near the critical point (t_e\approx t_c). In the thermodynamic limit (L\to \mathrm{\infty }), the bond- (or site-) addition dynamics can be viewed as a continuous realization of the growth of the correlation length \xi (t). The dynamic distribution {\mathcal{P}}_X(s) can therefore be regarded as a weighted superposition of the static distributions P_X(s,\xi ) over all possible correlation lengths \xi (t) sampled during the evolution. Since the dynamic process proceeds uniformly in time, the weighting factor W(\xi ) for a static distribution P_X(s,\xi ) can be directly characterized by the corresponding time interval \mathrm{d}t, leading to

(16)W(\xi )\mathrm{d}\xi =\mathrm{d}t.

At each time t with \xi (t)\gg 1, the system behaves as if it were in a quasi-critical state characterized by the correlation length \xi \sim |t-t_c{|}^{-\nu }. Differentiating \xi \sim |t-t_c{|}^{-\nu } yields

(17)\mathrm{d}\xi \sim |t-t_c{|}^{-\nu -1}\mathrm{d}t

(18)\sim {\xi }^{1+1/\nu }\mathrm{d}t.

Substituting into Eq. (16), we have

(19)W(\xi )\sim {\xi }^{-1-1/\nu }.

Recasting this expression in terms of the cutoff size s_{\xi }\sim {\xi }^{d_f} gives

(20)W(s_{\xi })\sim s_{\xi }^{-\sigma -1}.

This shows that the weighting factor W(s_{\xi }) itself obeys scaling with respect to the cutoff size s_{\xi } or the correlation length \xi.

The dynamic distribution {\mathcal{P}}_X(s) is therefore obtained as the integral of the static distribution P_X(s,s_{\xi }) weighted by W(s_{\xi }), over all relevant cutoff sizes s_{\xi }\ge s,

(21){\mathcal{P}}_X(s)\sim {{\displaystyle \int }}_s^{\mathrm{\infty }}W(s_{\xi })P_X(s,s_{\xi })\mathrm{d}{s}_{\xi }

(22)\sim s^{-{\tau }_X}{{\displaystyle \int }}_s^{\mathrm{\infty }}{s}_{\xi }^{-\sigma -1}{\tilde{P}}_X(s/s_{\xi })\mathrm{d}{s}_{\xi }.

This integral is simplified by the asymptotic behavior of the scaling function: {\tilde{P}}_X(s/s_{\xi }) approaches a constant near s_{\xi }=s and vanishes rapidly as s_{\xi }\to \mathrm{\infty }. Consequently, we have

(23){\mathcal{P}}_X(s)\sim s^{-{\tau }_X-\sigma }.

Therefore, the dynamic Fisher exponents satisfy

(24){\tau }_X^{\prime }={\tau }_X+\sigma ,

confirming Eqs. (14) and (15).

3.3.2 Dynamic contribution of the giant cluster

When the dynamic process proceeds into the supercritical regime (t_e>t_c), a macroscopic (giant) cluster emerges and continues to grow. As a result, the overall dynamic cluster-size distribution {\mathcal{P}}_S(s) inevitably acquires contributions from this giant cluster when the dynamic statistics are accumulated over the broader critical region (\xi \gg 1). In contrast, by construction, the gap distribution {\mathcal{P}}_G(s) naturally excludes the giant cluster and therefore retains the scaling form derived above.

In the thermodynamic limit (L\to \mathrm{\infty }), the number density of the giant cluster in a static configuration at t>t_c can be modeled by a Dirac delta function

(25)n_1(s)={\displaystyle \frac{1}{V}}\delta (s-C_1),

where C_1 is the size of the largest (giant) cluster, and V=L^d is the system volume. The normalization V\int n_1(s)\mathrm{d}s=1 ensures the uniqueness of the giant cluster. For a randomly inserted bond or site, the probability that it belongs to the giant cluster is m\equiv C_1/V, yielding the static size distribution of the giant cluster

(26)\begin{array}{rl}{P}_1(s)& =ms{n}_1(s)\\ & =\frac{s{C}_1}{V^2}\delta (s-C_1).\end{array}

Here, the additional factor s accounts for the fact that, when a bond or site is chosen at random, the probability of sampling a cluster of size s is proportional to its size – larger clusters contain more sites or bonds and are thus more likely to be selected. This is consistent with the general relation P_S(s)=sn(s) introduced earlier.

The dynamic process beyond t_c can be viewed as a temporal accumulation of configurations with an increasing giant cluster size C_1(t). The corresponding dynamic distribution {\mathcal{P}}_1(s) can then be obtained by integrating P_1(s) over C_1\in [1,V\to \mathrm{\infty }) with an C_1-dependent weight W(C_1). Similar to Eq. (16), W(C_1) is also characterized by the time interval \mathrm{d}t,

(27)W(C_1)\mathrm{d}{C}_1\sim \mathrm{d}t.

Note that the probability m(t)\equiv C_1(t)/V is just the order parameter of percolation transition. When \xi \gg 1, it obeys the scaling m(t)\sim (t-t_c{)}^{\beta } for t>t_c, so that C_1(t)\sim V(t-t_c{)}^{\beta }. This yields

(28)\mathrm{d}{C}_1\sim V(t-t_c{)}^{\beta -1}\mathrm{d}t

(29)\sim V^{1/\beta }{C}_1^{1-1/\beta }\mathrm{d}t,

leading to

(30)W(C_1)\sim V^{-1/\beta }{C}_1^{1/\beta -1}.

Therefore, the dynamic distribution {\mathcal{P}}_1(s) of the giant cluster can be expressed as

(31)\begin{array}{rl}{\mathcal{P}}_1(s)& \sim {\int }_1^{\mathrm{\infty }}W(C_1)P_1(s)\mathrm{d}{C}_1\\ & \sim V^{-2-1/\beta }s{\int }_1^{\mathrm{\infty }}{C}_1^{1/\beta }\delta (s-C_1)\mathrm{d}{C}_1\\ & \sim V^{-2-1/\beta }{s}^{1+1/\beta }.\end{array}

This expression applies only to the giant-cluster sector s>s_{\xi }, where s is of order C_1.

Consequently, the total dynamic cluster-size distribution under t_e>t_c can be expressed as a composite form

(32){\mathcal{P}}_S(s)=s^{-(\tau +\sigma -1)}{\tilde{\mathcal{P}}}_S(s/s_{\xi })+V^{-2-1/\beta }{s}^{1+1/\beta }{\tilde{\mathcal{P}}}_1(s/s_{\xi }),

where the scaling function {\tilde{\mathcal{P}}}_S(s/s_{\xi }) is nearly constant for and vanishes rapidly otherwise, while {\tilde{\mathcal{P}}}_1(s/s_{\xi }) behaves oppositely. Physically, the two contributions in Eq. (32) are statistically dominated by different stages of the dynamic process. The first term is dominated by configurations near and below t_c, since for t>t_c a randomly added bond or site is increasingly likely to involve the giant cluster, making events that connect only finite clusters rare. In contrast, the second term arises entirely from the supercritical regime and reflects the temporal accumulation of the growing giant cluster.

When s>s_{\xi }, the first term in Eq. (32) vanishes. As a result, for t_e>t_c, the tail of the dynamic cluster-size distribution is governed by the contribution of the giant cluster and behaves as Eq. (31). For a finite system with fixed volume V, this implies an algebraic increase {\mathcal{P}}_S(s)\sim s^{1+1/\beta } in the large-s regime. Consequently, {\mathcal{P}}_S(s) grows with s, giving rise to the hump observed at large s in Fig. 1(b1). For 2D percolation, \beta =5/36, yielding {\mathcal{P}}_S(s)\sim s^{41/5}, in agreement with the numerical results in Fig. 1(b1). This growth persists until s approaches the system-size cutoff s\sim V, corresponding to the largest possible cluster in the dynamic process. Substituting s\sim V into Eq. (31) yields {\mathcal{P}}_S(s)\sim V^{-1}\sim s^{-1}, indicating that the height of the hump decays as \sim s^{-1}. Notably, this decay is independent of the critical exponents and thus universal, as confirmed by Fig. 1(b1) and additional examples presented later.

4 Self-similar dynamics in broad percolation systems

Beyond providing an alternative route to reproduce known critical properties, the dynamic approach offers a general framework to characterize criticality from a temporal perspective. In particular, it is applicable to systems where the static ensemble description is problematic or where the critical point p_c cannot be precisely determined. In the following, we present representative examples that illustrate the generality and robustness of this dynamic approach.

4.1 Site percolation in 2D

For site percolation, the dynamic process proceeds by occupying sites randomly one by one. Unlike bond percolation, adding a single site can merge multiple clusters simultaneously. Following Eq. (4), the gap is defined as the total size of all merged clusters excluding the largest one.

Figure 2(a) shows simulation results on square lattices. Both the dynamic gap-size distribution, {\mathcal{P}}_G(s,L)\sim s^{-{\tau }_G^{\prime }}, and the dynamic cluster-size distribution, {\mathcal{P}}_S(s,L)\sim s^{-{\tau }_S^{\prime }} (inset), exhibit clear power-law scaling with exponents {\tau }_G^{\prime }=187/91 and {\tau }_S^{\prime }=132/91. These values coincide with those obtained for bond percolation in 2D (Fig. 1), confirming the universality of the dynamic scaling across site and bond percolation. The inset further shows that, for t_e=1, a bump appears for large s, contributed by the giant cluster. As predicted by Eq. (31), a scaling \sim s^{1+1/\beta }\sim s^{41/5} exists for large s, with the bump height scaling as \sim s^{-1}.

4.2 Bond percolation in 1D

In 1D percolation, the threshold is trivially located at t_c=1, where all sites become connected, so no conventional critical behavior is observed. Nevertheless, as shown in Fig. 2(b), the dynamically obtained gap-size distribution during the bond-insertion process still follows a clear power law, {\mathcal{P}}_G(s,L)\sim s^{-{\tau }_G^{\prime }}, with {\tau }_G^{\prime }=2, consistent with the known Fisher exponent \tau =1+d/d_f=2. This reflects the continuous nature of the transition with 1/\nu =d_f=1.

Since \sigma =1 in 1D, the corresponding dynamic cluster-size distribution {\mathcal{P}}_S(s,L) exhibits the same Fisher exponent, {\tau }_S^{\prime }=\tau -1+\sigma =2, as shown in the inset of Fig. 2(b). Because t_c=1, the system remains effectively subcritical even as t_e=1, and thus the bump in {\mathcal{P}}_S(s,L) at large s scales as \sim s^{-{\tau }_S^{\prime }}\sim s^{-2}, in contrast to the scaling \sim s^{-1} observed in models with .

4.3 Bond percolation on scale-free networks

A scale-free (SF) network is characterized by a degree distribution p_k\sim k^{-\lambda }, where k denotes the degree of a randomly chosen node, i.e., the number of links connected to the node. For small values of \lambda, the network is highly heterogeneous and contains hubs with very large degrees, whereas for large \lambda the degree distribution decays rapidly, and high-degree nodes become increasingly rare. Percolation on SF networks exhibits critical behavior that depends sensitively on the degree exponent \lambda [45–47]. However, it is difficult to achieve precise numerical confirmation, due to strong finite-size corrections [48, 49]. In simulations, SF networks can be generated using the configuration model [50, 51], in which each node is first assigned a degree drawn randomly from degree distribution p_k, and links are then randomly connected under the constraints of no self-loops or multiple edges.

As an example, we focus on \lambda >3, where a finite-threshold continuous transition exists and the critical exponents are well established. Analytically, \tau =(2\lambda -3)/(\lambda -2) and \sigma =(\lambda -3)/(\lambda -2) [45–47]; for \lambda =3.5, this gives \tau =8/3 and \sigma =1/3. As shown in Fig. 3(a), the dynamically measured gap-size distribution follows {\mathcal{P}}_G(s,L)\sim s^{-{\tau }_G^{\prime }} with {\tau }_G^{\prime }=8/3, in agreement with the theoretical Fisher exponent \tau.

For the dynamic cluster-size distribution {\mathcal{P}}_S(s,L)\sim s^{-{\tau }_S^{\prime }}, substituting the exponents \tau =(2\lambda -3)/(\lambda -2) and \sigma =(\lambda -3)/(\lambda -2) into {\tau }_S^{\prime }=\tau -1+\sigma yields {\tau }_S^{\prime }=2, independent of \lambda. The numeric result of the case \lambda =3.5 in the inset of Fig. 3(a) is well consistent with this scaling. It is pointed out that this exponent {\tau }_S^{\prime }=2 coincides with the mean-field result for percolation on complete graphs, for which the mean-field values \tau =5/2 and \sigma =1/2 also give {\tau }_S^{\prime }=2. This indicates that certain mean-field characteristics emerge for \lambda >3. Note that it is generally recognized that percolation on SF networks with \lambda \ge 4 exhibits the standard mean-field critical behavior [45].

Furthermore, from Eq. (31), the large-s regime in {\mathcal{P}}_S(s,L) follows the scaling \sim s^{1+1/\beta }, which reduces to \sim s^{\lambda -2} with \beta =1/(\lambda -3). The black dashed line in the inset of Fig. 3(b) clearly demonstrates this scaling \sim s^{1.5} for \lambda =3.5. The scaling \sim s^{-1} for the height of humps is also confirmed in the inset (blue dashed line).

4.4 Explosive percolation on complete graphs

We next apply the gap-dynamics framework to EP on complete graphs [8, 14–19, 40, 41]. EP follows the Achlioptas process [8], where at each step two candidate bonds are randomly selected and the bond minimizing the product of the sizes of the clusters it connects is added. Static FSS analyses often display anomalies reminiscent of discontinuous transitions; for a review, see Ref. [52]. Although EP is now understood to be continuous, extracting critical behaviors accurately remains challenging. Recently, two of us (ML and YD) have shown that EP obeys standard FSS when an event-based ensemble is applied [40, 41].

Applying the gap-dynamics framework, both the gap-size and cluster-size distributions exhibit clear power-law behavior [Fig. 3(b)]. Using previously reported values for the fractal dimension d_f=0.935 and the correlation-length exponent 1/\nu =0.740 [40], which are defined with respect to the system volume V, we obtain \tau =1+1/d_f\approx 2.07 and \sigma =1/(\nu d_f)\approx 0.79. As shown in Fig. 3(b) and its inset, the dynamically extracted exponents, {\tau }_G^{\prime }\approx 2.07 and {\tau }_S^{\prime }\approx 1.86, satisfy the relations {\tau }_G^{\prime }=\tau and {\tau }_S^{\prime }=\tau +\sigma -1. Moreover, from \beta =(1-d_f)\nu, the scaling for large s is \sim s^{1+1/\beta }\sim s^{12.38}. As shown in the inset of Fig. 3(b), there is good agreement between the simulation results and both the scaling \sim s^{12.38} at large s (black dashed line) and the decay \sim s^{-1} of the hump height (blue dashed line).

These results demonstrate the validity and robustness of the dynamic approach, even in systems exhibiting anomalous scaling under the conventional ensemble. Notably, the bond-insertion rule in EP endows it with an intrinsically dynamic character, making the dynamic ensemble a natural and informative perspective.

4.5 Rigidity percolation in 2D

Rigidity percolation follows the same bond-insertion rule as ordinary bond percolation, but focuses on rigid clusters capable of transmitting mechanical stress. It has been applied to glasses, gels, amorphous solids, jamming, fibrous networks, and living tissues [53–63]. Despite its broad relevance, critical properties of rigidity percolation remain elusive, with debated critical exponents and universality classes [10–13, 20–25].

In addition to exhibiting well-defined gap scaling as in the above percolation models, the gap-dynamics framework here serves mainly as an illustrative example; a detailed analysis of rigidity percolation is presented in Ref. [31]. In rigidity percolation, the gap-dynamics framework uncovers a cascade cluster-merging process: the insertion of a single bond can rigidify a nonrigid subgraph by activating multiple previously nonrigid links, triggering a collective transition. This process is quantified by \mathcal{K}, the number of clusters merged when a bond is inserted. As shown in Fig. 4, the \mathcal{K} distribution {\mathcal{P}}_K(s,L)\sim s^{-{\tau }_K} in an entire dynamic process follows a clear power law with {\tau }_K\approx 3.47, also revealing temporal self-similarity.

The average number of clusters merged, K(t)\equiv ⟨\mathcal{K}(t)⟩, exhibits a pronounced peak at the critical point t_c (inset of Fig. 4). While the peak does not diverge and curves for different L collapse without further rescaling, its location serves as a pseudocritical point. Analyzing the FSS of this pseudocritical point yields t_c=0.6602778(10), improving the precision by roughly two orders of magnitude compared to previous estimates t_c=0.6602(3) [10, 11].

5 Dynamic self-similarity in sandpile

The dynamic sampling method can also reveal nonequilibrium scaling behaviors in SOC systems. As a representative example, we consider the BTW sandpile model on a square lattice with open boundaries. Starting from an empty lattice, grains are added one by one to random sites. When a site accumulates four or more grains, it topples, redistributing grains to its four neighbors and possibly triggering further topplings. Grains can exit the system via the boundaries. After sufficiently long evolution, the system reaches a stationary critical state without any tuning parameter, characteristic of SOC.

Instead of focusing on the stationary regime as in conventional studies, we examine the initial nonequilibrium dynamics from the inception up to a pseudocritical time t_L, defined as the earliest avalanche that spans the system by reaching two opposite boundaries. At each time t\le t_L, we record the total number of topplings \mathcal{S} (toppled grains) and the toppling area \mathcal{A} (distinct toppled sites), and compute their dynamic distributions {\mathcal{P}}_S and {\mathcal{P}}_A. As shown in Figs. 5(a1) and (b1), both {\mathcal{P}}_S and {\mathcal{P}}_A exhibit power-law behavior \sim s^{-{\tau }^{\prime }}, sharing the same exponent {\tau }^{\prime }\approx 1.69. This value significantly differs from the stationary SOC exponents, which range from 1.05 to 1.29 [34, 35, 64–67], highlighting the distinct nature of the early-time dynamics.

We also analyze the avalanche diameter \mathcal{R}, defined as the span distance of toppled sites along one direction. The joint distribution \mathcal{P}(A,R) [inset of Fig. 5(b1)] reveals that \mathcal{A} scales as \sim {\mathcal{R}}^2, both in average and fluctuations, indicating a fractal dimension d_A=2. In contrast, the distribution \mathcal{P}(S,A) [inset of Fig. 5(a1)] broadens with increasing A: its lower boundary \mathcal{S}\sim \mathcal{A} corresponds to avalanches with single topplings per site, while the upper boundary \mathcal{S}\sim {\mathcal{A}}^{\kappa } with \kappa \gtrsim 1.34 captures cascades involving multiple topplings per site. If one instead considers only averaged quantities, S=⟨\mathcal{S}⟩ and A=⟨\mathcal{A}⟩, a simple scaling between S and A can still be observed, but it fails to reveal the correct underlying scaling behavior.

To further compare toppling area \mathcal{A} and toppling size \mathcal{S}, we examine their distributions P_A(s,L) and P_S(s,L) at t_L. For toppling area, plotting the rescaled variable x\equiv s/L^2 leads to excellent data collapse across different L [see Fig. 5(b2)], confirming that A is a well-defined fractal object with d_A=d. However, for toppling size, no single exponent d_S leads to full data collapse. As shown in Fig. 5(a2), to align the small-s regime, a small exponent d_S\approx 2 is needed, while a larger value d_S\approx 2.61 matches the peaks. Furthermore, the semi-log plot in Fig. 5(a2) suggests that the large-s tail decays exponentially as \sim {\mathrm{e}}^{-s/L^{d_S}}, indicating a non-power-law cutoff. Thus, S is not characterized by a single fractal dimension and cannot be viewed as a well-defined fractal object in the conventional sense.

In a conventional analysis for the SOC regime [36], Tebaldi et al. performed FSS of the avalanche-size moments ⟨{\mathcal{S}}^q⟩\sim L^{\theta (q)} and found \theta (1)=2, \theta (2)-\theta (1)=2.7, and \theta (9)-\theta (8)=2.9, suggesting that the size statistics exhibit weak multifractality. Our dynamic distributions \mathcal{P}(S,A) and {\mathcal{P}}_S(s,L) offer a more direct and vivid manifestation of this multifractal behavior. Based on the exponential cutoff form, we further conjecture that the maximal cutoff exponent for S is d_S=3.

6 Discussion

We reveal robust temporal self-similarity in various dynamic observables in a wide variety of percolation systems and BTW sandpile models. The self-similar dynamics not only captures scaling behaviors that were thought to emerge only in critical spatial structures, but also reveals rich critical phenomena that are challenging to be extracted by static analysis. Specifically, we theoretically established the dynamic scaling relations, {\tau }_G^{\prime }=\tau and {\tau }_S^{\prime }=\tau +\sigma -1, linking the dynamic Fisher exponents directly to the static exponents.

The dynamic method serves as a powerful new paradigm for extracting critical behavior by incorporating both temporal evolution and spatial information. First, dynamic self-similarity integrates critical properties over a broad parameter range, thereby avoiding the need for precise identification of the critical point. Second, beyond temporal observables, spatial quantities sampled at the pseudocritical point t_L can be analyzed using FSS theory, and the flexibility in defining t_L provides an additional degree of freedom to tailor the analysis and mitigate finite-size corrections. Important applications include high-dimensional percolation [42], long-ranged percolation [68], and percolation on SF networks [49]. Third, dynamic approaches can reveal critical behavior that remains elusive under static averaging. A representative example is EP, whose critical behavior has been clarified by dynamic analyses [40, 41, 69]. Fourth, dynamic analysis may encode richer information than static snapshots, such as the cascade effect in rigidity percolation, and the weak multifractal character in the BTW sandpile model.

More broadly, the dynamic formulation provides a natural framework for systems whose critical behavior emerges through continuous evolution rather than well-defined steady states. It is readily applicable to a wide range of models with complex dynamics, including jamming and cascading failures, and offers promising perspectives for studying infrastructure resilience, biological systems, and geophysical processes [70], where static descriptions are often insufficient.

References

[1] S. K. Ma, Modern Theory of Critical Phenomena, Routledge, 2018

[2] D. Stauffer,A. Aharony, Introduction to Percolation Theory, 2nd Ed., Taylor & Francis, London, 1991

[3] B. Nienhuis , E. K. Riedel , and M. Schick , Magnetic exponents of the two-dimensional q-state Potts model, J. Phys. Math. Gen. 13(6), L189 (1980)

[4] B. Nienhuis , Critical behavior of two-dimensional spin models and charge asymmetry in the Coulomb gas, J. Stat. Phys. 34(5−6), 731 (1984)

[5] B. Nienhuis, Coulomb gas formulation of 2D phase transition, in: Phase transition and critical phenomena, edited by C. Domb and J. L. Lebowitz, 1987, Vol. 11

[6] J. L. Cardy, Conformal invariance, in: Phase Transition and Critical Phenomena, edited by C. Domb, Vol. 11, p. 55, Academic Press, London, 1987

[7] S. Smirnov and W. Werner , Critical exponents for two-dimensional percolation, Math. Res. Lett. 8(6), 729 (2001)

[8] D. Achlioptas , R. M. D’Souza , and J. Spencer , Explosive percolation in random networks, Science 323(5920), 1453 (2009)

[9] O. Riordan and L. Warnke , Explosive percolation is continuous, Science 333(6040), 322 (2011)

[10] D. J. Jacobs and M. F. Thorpe , Generic rigidity percolation: The pebble game, Phys. Rev. Lett. 75(22), 4051 (1995)

[11] D. J. Jacobs and M. F. Thorpe , Generic rigidity percolation in two dimensions, Phys. Rev. E 53(4), 3682 (1996)

[12] C. Moukarzel and P. M. Duxbury , Comparison of rigidity and connectivity percolation in two dimensions, Phys. Rev. E 59(3), 2614 (1999)

[13] M. A. Brière , M. V. Chubynsky , and N. Mousseau , Self-organized criticality in the intermediate phase of rigidity percolation, Phys. Rev. E 75(5), 056108 (2007)

[14] E. J. Friedman and A. S. Landsberg , Construction and analysis of random networks with explosive percolation, Phys. Rev. Lett. 103(25), 255701 (2009)

[15] R. M. Ziff, Explosive growth in biased dynamic percolation on two-dimensional regular lattice networks, Phys. Rev. Lett. 103(4), 045701 (2009)

[16] R. M. D’Souza and M. Mitzenmacher , Local cluster aggregation models of explosive percolation, Phys. Rev. Lett. 104(19), 195702 (2010)

[17] F. Radicchi and S. Fortunato , Explosive percolation: A numerical analysis, Phys. Rev. E 81(3), 036110 (2010)

[18] P. Grassberger , C. Christensen , G. Bizhani , S. W. Son , and M. Paczuski , Explosive percolation is continuous, but with unusual finite size behavior, Phys. Rev. Lett. 106(22), 225701 (2011)

[19] R. M. D’Souza and J. Nagler , Anomalous critical and supercritical phenomena in explosive percolation, Nat. Phys. 11(7), 531 (2015)

[20] S. Feng and P. N. Sen , Percolation on elastic networks: New exponent and threshold, Phys. Rev. Lett. 52(3), 216 (1984)

[21] Y. Kantor and I. Webman , Elastic properties of random percolating systems, Phys. Rev. Lett. 52(21), 1891 (1984)

[22] A. Hansen and S. Roux , Universality class of central-force percolation, Phys. Rev. B 40(1), 749 (1989)

[23] M. Plischke , Rigidity of disordered networks with bond-bending forces, Phys. Rev. E 76(2), 021401 (2007)

[24] L. Zhang , D. Z. Rocklin , B. G. Chen , and X. Mao , Rigidity percolation by next-nearest-neighbor bonds on generic and regular isostatic lattices, Phys. Rev. E. 91(3), 032124 (2015)

[25] W. G. Ellenbroek , V. F. Hagh , A. Kumar , M. Thorpe , and M. van Hecke , Rigidity loss in disordered systems: Three scenarios, Phys. Rev. Lett. 114(13), 135501 (2015)

[26] M. Lu,Y. F. Song,M. Li,Y. Deng, Self-similar gap dynamics in percolation and rigidity percolation, arXiv: 2024)

[27] S. V. Buldyrev , R. Parshani , G. Paul , H. E. Stanley , and S. Havlin , Catastrophic cascade of failures in interdependent networks, Nature 464(7291), 1025 (2010)

[28] J. Fan , J. Meng , Y. Ashkenazy , S. Havlin , and H. J. Schellnhuber , Climate network percolation reveals the expansion and weakening of the tropical component under global warming, Proc. Natl. Acad. Sci. USA 115(52), E12128 (2018)

[29] D. Li , B. Fu , Y. Wang , G. Lu , Y. Berezin , H. E. Stanley , and S. Havlin , Percolation transition in dynamical traffic network with evolving critical bottlenecks, Proc. Natl. Acad. Sci. USA 112(3), 669 (2015)

[30] J. Xie , F. Meng , J. Sun , X. Ma , G. Yan , and Y. Hu , Detecting and modelling real percolation and phase transitions of information on social media, Nat. Hum. Behav. 5(9), 1161 (2021)

[31] M. Lu,Y. F. Song,Q. Shi,M. Li,Y. Deng, High-precision dynamic Monte Carlo study of rigidity percolation, arXiv: 2026)

[32] P. Bak , C. Tang , and K. Wiesenfeld , Self-organized criticality: An explanation of the 1/f noise, Phys. Rev. Lett. 59(4), 381 (1987)

[33] C. Tang and P. Bak , Critical exponents and scaling relations for self-organized critical phenomena, Phys. Rev. Lett. 60(23), 2347 (1988)

[34] M. Engsig and K. Sneppen , Fractals in the critical attractor of the classical sandpile model, Phys. Rev. Lett. 134(18), 187201 (2025)

[35] S. S. Manna , Describing self-organized criticality as a continuous phase transition, Phys. Rev. E 111(2), 024111 (2025)

[36] C. Tebaldi , M. De Menech , and A. L. Stella , Multifractal scaling in the Bak–Tang–Wiesenfeld sandpile and edge events, Phys. Rev. Lett. 83(19), 3952 (1999)

[37] M. E. J. Newman and R. M. Ziff , Fast Monte Carlo algorithm for site or bond percolation, Phys. Rev. E 64(1), 016706 (2001)

[38] S. N. Dorogovtsev , J. F. F. Mendes , and A. N. Samukhin , Anomalous percolation properties of growing networks, Phys. Rev. E 64(6), 066110 (2001)

[39] B. Bollobas , S. Janson , and O. Riordan , The phase transition in the uniformly grown random graph has infinite order, Random Structures Algorithms 26(1−2), 1 (2004)

[40] M. Li , J. Wang , and Y. Deng , Explosive percolation obeys standard finite-size scaling in an event-based ensemble, Phys. Rev. Lett. 130(14), 147101 (2023)

[41] M. Li , J. Wang , and Y. Deng , Explosive percolation in finite dimensions, Phys. Rev. Res. 6(3), 033319 (2024)

[42] M. Li,S. Fang,J. Fan,Y. Deng, Crossover finite-size scaling theory and its applications in percolation, arXiv: 2024)

[43] J. Fan , J. Meng , Y. Liu , A. A. Saberi , J. Kurths , and J. Nagler , Universal gap scaling in percolation, Nat. Phys. 16(4), 455 (2020)

[44] Y. Deng , W. Zhang , T. M. Garoni , A. D. Sokal , and A. Sportiello , Some geometric critical exponents for percolation and the random-cluster model, Phys. Rev. E 81(2), 020102 (2010)

[45] R. Cohen , D. ben-Avraham , and S. Havlin , Percolation critical exponents in scale-free networks, Phys. Rev. E 66(3), 036113 (2002)

[46] D. S. Lee , K. I. Goh , B. Kahng , and D. Kim , Evolution of scale-free random graphs: Potts model formulation, Nucl. Phys. B 696(3), 351 (2004)

[47] L. Cirigliano , G. Timar , and C. Castellano , Scaling and universality for percolation in random networks: A unified view, Phys. Rev. E 110(6), 064303 (2024)

[48] Z. Wu , C. Lagorio , L. A. Braunstein , R. Cohen , S. Havlin , and H. E. Stanley , Numerical evaluation of the upper critical dimension of percolation in scale-free networks, Phys. Rev. E 75(6), 066110 (2007)

[49] X. Zhao , L. Yang , D. Peng , R. R. Liu , and M. Li , Finite-size scaling of percolation on scale-free networks, Chaos Solitons Fractals 200, 117076 (2025)

[50] M. Molloy,B. Reed, A critical point for random graphs with a given degree sequence, Random Structures Algorithms 6(2–3), 161 (1995)

[51] M. E. J. Newman , S. H. Strogatz , and D. J. Watts , Random graphs with arbitrary degree distributions and their applications, Phys. Rev. E 64(2), 026118 (2001)

[52] R. M. D’Souza , J. Gomez-Gardees , J. Nagler , and A. Arenas , Explosive phenomena in complex networks, Adv. Phys. 68(3), 123 (2019)

[53] M. F. Thorpe,D. J. Jacobs,M. V. Chubynsky,J. C. Phillips, Self-organization in network glasses, J. Non-Cryst. Solids 266–269, 859 (2000)

[54] C. P. Broedersz , X. Mao , T. C. Lubensky , and F. C. MacKintosh , Criticality and isostaticity in fibre networks, Nat. Phys. 7(12), 983 (2011)

[55] S. Zhang , L. Zhang , M. Bouzid , D. Z. Rocklin , E. Del Gado , and X. Mao , Correlated rigidity percolation and colloidal gels, Phys. Rev. Lett. 123(5), 058001 (2019)

[56] J. Rouwhorst , C. Ness , S. Stoyanov , A. Zaccone , and P. Schall , Nonequilibrium continuous phase transition in colloidal gelation with short-range attraction, Nat. Commun. 11(1), 3558 (2020)

[57] J. Rouwhorst , P. Schall , C. Ness , T. Blijdenstein , and A. Zaccone , Nonequilibrium master kinetic equation modeling of colloidal gelation, Phys. Rev. E 102(2), 022602 (2020)

[58] A. Zaccone and E. Scossa-Romano , Approximate analytical description of the nonaffine response of amorphous solids, Phys. Rev. B. 83(18), 184205 (2011)

[59] S. Henkes , D. A. Quint , Y. Fily , and J. Schwarz , Rigid cluster decomposition reveals criticality in frictional jamming, Phys. Rev. Lett. 116(2), 028301 (2016)

[60] H. Dashti , A. A. Saberi , S. Rahbari , and J. Kurths , Emergence of rigidity percolation in flowing granular systems, Sci. Adv. 9(35), eadh5586 (2023)

[61] H. A. Vinutha , F. D. Diaz Ruiz , X. Mao , B. Chakraborty , and E. Del Gado , Stress–stress correlations reveal force chains in gels, J. Chem. Phys. 158(11), 114104 (2023)

[62] N. I. Petridou , B. Corominas-Murtra , C. P. Heisenberg , and E. Hannezo , Rigidity percolation uncovers a structural basis for embryonic tissue phase transitions, Cell 184(7), 1914 (2021)

[63] J. Rozman , M. Krajnc , and P. Ziherl , Basolateral mechanics prevents rigidity transition in epithelial monolayers, Phys. Rev. Lett. 133(16), 168401 (2024)

[64] S. Lübeck and K. D. Usadel , Numerical determination of the avalanche exponents of the Bak–Tang–Wiesenfeld model, Phys. Rev. E 55(4), 4095 (1997)

[65] L. P. Kadanoff , S. R. Nagel , L. Wu , and S. Zhou , Scaling and universality in avalanches, Phys. Rev. A 39(12), 6524 (1989)

[66] S. Manna , Critical exponents of the sandpile models in two dimensions, Physica A 179(2), 249 (1991)

[67] E. Milshtein , O. Biham , S. Solomon , and Universality classes in isotropic , Abelian, and non-Abelian sandpile models, Phys. Rev. E 58(1), 303 (1998)

[68] Z. Liu,T. Xiao,Z. Fan,Y. Deng, Two-dimensional percolation model with long-range interaction, arXiv: 2025)

[69] L. Yang and M. Li , Emergence of biconnected clusters in explosive percolation, Phys. Rev. E 110(1), 014122 (2024)

[70] S. Fang,Q. Lin,J. Meng,B. Chen,J. Nagler,Y. Deng,J. Fan, Universal scaling of gap dynamics in percolation, arXiv: 2024)

RIGHTS & PERMISSIONS

Higher Education Press