Newman-Ziff algorithm for the bootstrap percolation: application to the Archimedean lattices

TL;DR

This paper introduces efficient algorithms based on the Newman-Ziff method to analyze bootstrap and diffusion percolation models on Archimedean lattices, revealing they share the same universality class as classical percolation.

Contribution

The authors extend the Newman-Ziff algorithm to bootstrap and diffusion percolation, enabling high-precision calculation of thresholds and critical exponents on Archimedean lattices.

Findings

Bootstrap and diffusion percolation have the same critical exponents as classical percolation.

The algorithms allow precise determination of percolation thresholds.

Percolation transitions are continuous in these models.

Abstract

We propose very efficient algorithms for the bootstrap percolation and the diffusion percolation models by extending the Newman-Ziff algorithm of the classical percolation [M. E. J. Newman and R. M. Ziff, Phys. Rev. Lett. 85 (2000) 4104]. Using these algorithms and the finite-size-scaling, we calculated with high precision the percolation threshold and critical exponents in the eleven two-dimensional Archimedean lattices. We present the condition for the continuous percolation transition in the bootstrap percolation and the diffusion percolation, and show that they have the same critical exponents as the classical percolation within error bars in two dimensions. We conclude that the bootstrap percolation and the diffusion percolation almost certainly belong to the same universality class as the classical percolation.

Tables

Table 1. Table 1: Name, coordination number ( Δ Δ \Delta ), percolation threshold ( p BP m superscript subscript 𝑝 BP 𝑚 p_{\mathrm{\scriptscriptstyle BP}}^{m} and p DP k superscript subscript 𝑝 DP 𝑘 p_{\mathrm{\scriptscriptstyle DP}}^{k} ), and critical exponents ( ν 𝜈 \nu and β 𝛽 \beta ) for the BP and DP with the continuous percolation transition of the eleven Archimedean lattices.

		Bootstrap percolation				Diffusion percolation
Lattice	$\Delta$	$m$	$p_{\mathrm{BP}}^m$	$\nu$	$\beta /\nu$	$k$	$p_{\mathrm{DP}}^k$	$\nu$	$\beta /\nu$
T1	6	$1$	0.49997(4)	1.336(3)	0.104(1)	$6$	0.50000(1)	1.335(1)	0.104(1)
($3^6$)		$2$	0.49999(4)	1.336(3)	0.104(1)	$5$	0.49999(2)	1.334(1)	0.104(1)
		$3$	0.62915(5)	1.335(3)	0.104(2)	$4$	0.37083(4)	1.334(2)	0.103(1)
T2	4	$1$	0.59272(3)	1.335(4)	0.104(1)	$4$	0.54731(1)	1.333(4)	0.104(1)
($4^4$)		$2$	0.59272(3)	1.336(3)	0.104(1)	$3$	0.42037(2)	1.333(2)	0.104(1)
T3	3	$1$	0.69710(5)	1.337(3)	0.104(1)	$3$	0.56008(1)	1.333(1)	0.104(1)
($6^3$)		$2$	0.69703(4)	1.335(3)	0.104(1)	$2$	0.30943(2)	1.333(2)	0.104(1)
T4	5	$1$	0.57948(3)	1.337(3)	0.104(1)	$5$	0.57950(1)	1.332(2)	0.104(1)
($3^4,6$)		$2$	0.57948(3)	1.336(3)	0.104(1)	$4$	0.48450(1)	1.334(1)	0.104(1)
		$3$	0.73227(3)	1.333(2)	0.103(3)	$3$	0.26936(3)	1.334(1)	0.104(1)
T5	5	$1$	0.55020(3)	1.337(3)	0.104(1)	$5$	0.54387(1)	1.333(1)	0.104(1)
($3^3,4^2$)		$2$	0.55020(4)	1.335(3)	0.104(1)	$4$	0.47648(2)	1.333(2)	0.104(1)
		$3$	0.71884(4)	1.330(4)	0.103(1)	$3$	0.28165(1)	1.332(1)	0.103(1)
T6	5	$1$	0.55080(3)	1.337(4)	0.104(1)	$5$	0.54108(1)	1.333(1)	0.104(1)
($3^2,4,3,4$)		$2$	0.55080(4)	1.336(3)	0.104(1)	$4$	0.47072(1)	1.334(1)	0.104(1)
		$3$	0.72813(5)	1.336(4)	0.101(4)	$3$	0.27194(1)	1.335(1)	0.103(1)
T7	4	$1$	0.62180(3)	1.337(3)	0.104(1)	$4$	0.57502(1)	1.332(2)	0.104(1)
($3,4,6,4$)		$2$	0.62178(4)	1.336(3)	0.104(1)	$3$	0.42652(1)	1.333(1)	0.104(1)
		$3$	0.86713(5)	1.337(4)	0.103(3)	$2$	0.13447(3)	1.334(2)	0.104(1)
T8	4	$1$	0.65268(3)	1.336(5)	0.104(1)	$4$	0.58661(1)	1.334(2)	0.104(1)
($3,6,3,6$)		$2$	0.65269(4)	1.337(4)	0.104(1)	$3$	0.39451(2)	1.335(2)	0.104(1)
T9	3	$1$	0.80787(3)	1.338(4)	0.104(1)	$3$	0.65335(1)	1.332(1)	0.104(1)
($3,{12}^2$)		$2$	0.80787(3)	1.338(3)	0.104(1)	$2$	0.34028(4)	1.333(2)	0.104(1)
T10	3	$1$	0.74779(3)	1.339(4)	0.104(1)	$3$	0.61644(1)	1.332(1)	0.104(1)
($4,6,12$)		$2$	0.74779(3)	1.337(4)	0.104(1)	$2$	0.31816(3)	1.333(1)	0.104(1)
T11	3	$1$	0.72971(3)	1.336(5)	0.104(1)	$3$	0.58862(1)	1.334(2)	0.104(1)
($4,8^2$)		$2$	0.72971(4)	1.336(5)	0.104(1)	$2$	0.30280(4)	1.334(1)	0.104(1)

Equations

$$⟨Q(p)⟩=\sum _{n=0}^N\frac{N!}{n!(N-n)!}{p}^n(1-p{)}^{N-n}⟨Q(n)⟩,$$

$$\frac{d⟨Q(p)⟩}{dp}=\frac{d}{dp}\left[\sum _{n=0}^N\frac{N!}{n!(N-n)!}{p}^n(1-p{)}^{N-n}⟨Q(n)⟩\right]$$

$$=\sum _{n=0}^N\frac{N!}{n!(N-n)!}{p}^{n-1}(1-p{)}^{N-n-1}(n-Np)⟨Q(n)⟩$$

$$\mathrm{Max}\left[d{P}_{w1}(p,L)/dp\right]\sim L^{1/\nu },$$

$$\mathrm{Max}\left[d{P}_{w2}(p,L)/dp\right]\sim L^{1/\nu },$$

$${\left[d{P}_{w1}(p,L)/dp\right]}_{p=p_c^{(\infty )}}\sim L^{1/\nu },$$

$${\left[d{P}_{w2}(p,L)/dp\right]}_{p=p_c^{(\infty )}}\sim L^{1/\nu },$$

$$\text{mbox}and{P}_{\infty }(p_c^{(\infty )},L)\sim L^{-\beta /\nu }.$$

Full text

1 Introduction

The bootstrap percolation (BP) model [1] has attracted continuous attention for its various applications such as disordered dilute magnetic systems [1, 2], neuronal activity [3], jamming transition [4], and diffusion of innovations [5]. The BP process operates as follows: (i) Each site is occupied with the probability $p$ and empty otherwise; (ii) Occupied sites that have less than $m$ occupied neighbors become empty, and the process is repeated until all the occupied sites have at least $m$ occupied neighbors. The diffusion percolation (DP) model [6] is closely-related to the BP, and sometimes it is also called the BP [7, 8, 9]. In the DP, process (i) is the same as the BP, but process (ii) is different: empty sites that have at least $k$ occupied neighbors become occupied recursively, until all the empty sites have less than $k$ occupied neighbors.

The BP and DP have been studied on lattices [10, 11, 12, 13, 6, 14, 15, 16, 17, 18, 8], trees [1, 19, 20], and complex networks [21, 7, 22, 9]. A few facts are known about the BP and DP. Clearly, the BP with $m=0$ and DP with $k>\Delta_{\mathrm{max}}$ are the same as the classical percolation (CP) model, where $\Delta_{\mathrm{max}}$ is the maximum value of degree. For a given graph, the percolation threshold of the DP ($p_{\mathrm{\scriptscriptstyle DP}}^{k}$) is not larger than that of the CP ($p_{\mathrm{\scriptscriptstyle CP}}$), and that of the BP ($p_{\mathrm{\scriptscriptstyle BP}}^{m}$) is not less than that of the CP. The BP with $m=1$ or $m=2$ has the same percolation threshold as the CP. As for $\Delta$-regular lattices, if $m+k=\Delta+1$, then there is a close relationship between the BP and DP: (1) The sum of percolation thresholds of the BP and its corresponding DP is the same or larger than 1 ($p_{\mathrm{\scriptscriptstyle BP}}^{m}+p_{\mathrm{\scriptscriptstyle DP}}^{\Delta+1-m}\geq 1$) [6] (The sum is 1 for self-matching lattices such as the triangular lattice [23].); (2) There are lattice-dependent parameters $m_{c}$ and $k_{c}$, where the BP of $m$ and DP of $k$ have first-order percolation transitions with the percolation thresholds $p_{\mathrm{\scriptscriptstyle BP}}^{m}=1$ and $p_{\mathrm{\scriptscriptstyle DP}}^{k}=0$, respectively, if and only if $m>m_{c}$ and $k<k_{c}$; the two parameters satisfy $m_{c}+k_{c}=\Delta+1$. For the square and triangular lattices, $m_{c}=\Delta/2$ [10, 24, 6].

About the critical exponents and the universality class of continuous percolation transitions of the BP and DP, there have been long debates. There are several critical exponents for the percolation, and in fact any two of them determine all of the rest exponents through scaling relations [25]. For the CP, they are universal and depend only on the spatial dimension. While the BP of $m=1$ is known to have the same critical exponents as the CP [14], it is not clear for the DP and the BP of $m\geq 2$. Kogut and Leath argued that $\beta$ depends on $m$ for the BP from Monte-Carlo simulations on the square, triangular, and cubic lattices [10], and renormalization group studies confirmed it [11, 13]. Adler conjectured that $\nu$ is universal but $\beta$ is non-universal in the BP and DP [6, 14, 15]. To the contrary, other studies, which include most recent simulations, insist that both $\nu$ and $\beta$ are universal [12, 16, 17, 18]. They argued that non-universality of previous works can be from the small size of clusters used in the simulations. However, we judge that the universality class of the BP and DP is not definitely clear yet. In Ref. [17], for example, the Fisher exponent $\tau$ for the DP with $k=4$ was calculated on the triangular lattice to be $\tau=2.03\pm 0.04$, which is consistent with $\tau=187/91\approx 2.055$ for the two-dimensional CP. At first glance, it looks like a good evidence to support the conclusion of the same universality; however, $\beta$ obtained by the scaling relation $\beta/\nu=2(\tau-2)/(\tau-1)$ in two dimensions is $\beta/\nu=0.06\pm 0.08$, which is consistent with $\beta/\nu=5/48\approx 0.104$ of the CP but the uncertainty is too large to make any conclusion. In addition, the BP and DP have been studied only in four kinds of lattices (the square, triangular, honeycomb, and cubic lattices).

In this paper, we introduce efficient algorithms for the BP and DP models and present much more precise results on eleven two-dimensional Archimedean lattices. We calculated percolation threshold and critical exponents ($\nu$ and $\beta$) to present positive evidences for that the BP and DP belong to the same universality class as the CP in two dimensions.

2 Methods

The classical site percolation model can be simulated simply by the following two steps: (i) Fill each site with independent probability $p$ and leave it empty with probability $1-p$; (ii) Identify all the connected clusters to check whether a percolating cluster exists. The simulation should be repeated many times to make a statistical average for a given $p$. This algorithm is simple but inefficient. In order to get results with different probability $p^{\prime}$, the whole simulation should be done again. It is more difficult to calculate the derivative of a quantity with respect to the probability $p$, because numerical differentiation inevitably gives large error [26]. A more efficient algorithm was proposed by Newman and Ziff [27, 28], and it has become a standard method in classical percolation studies. The Newman-Ziff algorithm consists of four steps: (i) Initially, all sites are empty; (ii) Choose an empty site randomly and fill it; (iii) Update the information of connected clusters to check whether percolation occurs; (iv) Repeat steps (ii) and (iii) until all the sites are occupied. An efficient algorithm to update the information of connected clusters (tree-based union/find algorithms) is also presented in Ref. [28]. An average $\langle Q(n)\rangle$ is obtained by repeating the whole steps, where $Q(n)$ is any quantity (e.g., size of the largest cluster) for a fixed number of occupied sites $n$. In one run of the Newman-Ziff algorithm, $Q(n_{1})$ is correlated with $Q(n_{2})$ inevitably, but values of $Q(n)$ of different runs are absolutely independent and so the statistical averaging of $\langle Q(n)\rangle$ has no problem. A value $\langle Q(p)\rangle$ for a fixed occupation probability $p$ can be obtained by the transformation of

[TABLE]

where $N$ is the total number of sites. Therefore, once $\langle Q(n)\rangle$ is obtained, $\langle Q(p)\rangle$ can be calculated for all values of $p$. Another advantage of the Newman-Ziff algorithm is that the derivative can be obtained through

[TABLE]

without numerical differentiation [29].

However, the Newman-Ziff algorithm cannot be used directly in the BP or DP models, because filling of each site depends on the local environment. As for the DP, the Newman-Ziff algorithm can be modified as follows.

(1) Initially, all sites are empty.

(2) Make an array of all the sites in random order.

(3) Get one site by the array. If the site is empty, fill it and other sites that have at least $k$ occupied neighbors, recursively. If the site is already occupied, do nothing.

(4) Update the information about connected clusters to check whether percolation occurs.

(5) Repeat steps (3) and (4) until all the sites are occupied.

The whole steps are repeated to make an average quantity $\langle Q(n)\rangle$, and the transformation of Eq. (1) gives $\langle Q(p)\rangle$. This algorithm is equivalent to the original DP with fixed $p$ because the final state does not depend on the sequence of the filling process once an initial occupation is determined. The BP model can be simulated by the same way as DP: Initially all sites are filled and sites are emptied one-by-one in random order. This algorithm is called the avalanching bootstrap percolation of the second kind (ABP2) [30]. However, the algorithm is inefficient because it is difficult to identify and update the cluster information during simulation. Therefore, we propose a more efficient algorithm for the BP model by introducing preoccupied state in addition to empty and occupied states. When a site is chosen in the Newman-Ziff algorithm, it is occupied only when there are at least $m$ occupied neighbors; otherwise, it is assigned to be preoccupied. Preoccupied state means that the site will be occupied after the condition is satisfied. The algorithm can be presented as follows.

(1) Initially, all sites are empty.

(2) Make an array of all the sites in random order.

(3-a) Get one site by the array. Set the site into the preoccupied state.

(3-b) Identify the connected cluster of preoccupied sites that includes the site. (To accelerate the simulation, preoccupied sites with less than $m$ occupied or preoccupied neighbors can be excluded from the cluster, because they are impossible to be filled.)

(3-c) Fill tentatively all the sites of the cluster.

(3-d) Within the cluster, set sites that have less than $m$ occupied neighbors into the preoccupied state again, recursively, until all the occupied sites have at least $m$ occupied neighbors.

(4) Update the information about connected clusters to check whether percolation occurs.

(5) Repeat steps (3) and (4) until all the sites are occupied.

All steps of this algorithm are identical with the Newman-Ziff algorithm of the CP except for step (3), and the efficient routine that updates the information of connected clusters can be used without modification. Although we focus on two-dimensional lattices in this paper, both of the algorithms for the DP and BP can also be applied to any dimensional systems and to complex networks.

The CPU time needed in these algorithms is proportional to the lattice size by $T_{\mathrm{cpu}}\sim N^{\lambda}$ with $\lambda\approx 1.2$. CPU time for the DP, CP, and BP ($m<3$) are of the same order of magnitude for the same value of $N$; in the case of the BP with $m=3$, CPU time requirement is about 10 times more than that of the CP. Three-dimensional lattices show the same behavior, but more calculations would be needed than two dimensions because the number of sites increases as $N\sim L^{3}$, where $L$ is linear size. For a lattice of $N=10^{6}$, it takes about one second of CPU time for one sweep except for the BP of $m=3$ by Intel(R) Xeon(R) CPU of 2.2 GHz. Note that the running time can be easily reduced by parallelizing the algorithm. In the case of the traditional brute force approach, which calculates $\langle Q(p)\rangle$ directly, the CPU time requirement depends on the number of occupied sites. On average, it takes about half of the Newman-Ziff algorithm because some part of the lattice is not filled at all. However, it gives only $\langle Q(p)\rangle$ at a specific value of $p$, and another independent calculation is needed to get $\langle Q(p^{\prime})\rangle$ for $p^{\prime}\neq p$. In addition, it is practically impossible to get its derivative $d\langle Q(p)\rangle/dp$ with high precision. Therefore, our new algorithms are much more efficient except when only $\langle Q(p)\rangle$ is to be calculated at a known specific $p$.

Using these algorithms, we calculated the strength of the largest cluster (the probability that a site belongs to the largest cluster; $P_{\infty}$), average cluster size excluding the largest one ($M^{\prime}_{1}$), percolation probability in any direction ($P_{w1}$) and in both directions ($P_{w2}$), and proportion of occupied sites ($P_{o}$) as a function of $p$ for the CP, BP, and DP.

In this work, we consider eleven Archimedean lattices, which are vertex-transitive graphs made in two dimensions by edge-to-edge tiling of regular polygons whose vertices are surrounded by the same sequence of polygons. There are only eleven Archimedean lattices [31] and they are typically used for systematic studies [32, 33, 34]. They are shown in Fig. 1. The periodic boundary condition is used and the percolation is defined by the existence of a cluster that wraps all the way around the lattice. The number of lattice sites studied in this work is from $N=1296$ to $N=2\,560\,000$.

The percolation thresholds ($p_{\mathrm{\scriptscriptstyle DP}}^{k}$ and $p_{\mathrm{\scriptscriptstyle BP}}^{m}$) of infinite lattices are determined by the finite-size-scaling. The percolation threshold estimate of a finite lattice is determined by the probabilities of initial filling ($p$) that give the maximums of physical quantities that show critical behavior or their derivatives. The percolation threshold can also be found by the crossing points of percolation probabilities of different lattice size [28]. We averaged the percolation threshold values obtained by these methods to get the final estimate. Critical exponents $\nu$ and $\beta$ are obtained by the derivative of percolation probability and $P_{\infty}$, respectively [35, 29]. The correction-to-scaling [36, 37] is ignored.

Most of the results were produced by using Mersenne twister pseudo-random-number generator (MT19937) [38], which was confirmed reliable in a site percolation problem [39]. We also confirmed that other pseudo-random-number generators give equivalent results.

3 Results

Figure 2 shows the strength of the largest cluster $P_{\infty}(p,L)$, average cluster size excluding the largest one $M^{\prime}_{1}(p,L)$, percolation probability in any direction $P_{w1}(p,L)$ and in both directions $P_{w2}(p,L)$, and proportion of occupied sites $P_{o}(p,L)$ for the DP ($k=4$), CP, and BP ($m=3$) in the triangular lattice. Parameter $L$ is the linear size of the lattice. They all show continuous phase transition, which becomes sharper as the lattice size increases. The proportion of occupied sites $P_{o}(p,L)$ is independent of lattice size and does not show any critical behavior at the percolation threshold. Equivalent behavior is observed in the BP and DP of the other values of $m\leq m_{c}$ and $k\geq k_{c}$, and in the other Archimedean lattices. The percolation threshold of infinite lattices [$p_{c}^{(\infty)}=\lim_{L\rightarrow\infty}p_{c}(L)$] is determined by the finite-size-scaling: $[p_{c}(L)-p_{c}^{(\infty)}]\sim L^{-a}$. The percolation threshold estimate for a finite lattice $p_{c}(L)$ is determined by the probabilities of initial filling ($p$) that give the maximum values of $dP_{\infty}(p,L)/dp$, $M^{\prime}_{1}(p,L)$, $[P_{w2}(p,L)-P_{w1}(p,L)]$, $dP_{w1}(p,L)/dp$, and $dP_{w2}(p,L)/dp$. The fitting parameter $a$ depend on lattice structure, the physical quantity measured, and percolation type [40]. Figure 3 confirms the scaling behavior very well. In the case of the BP of $m=3$, however, deviation from the scaling is large in small lattices for $dP_{\infty}(p,L)/dp$ and $M^{\prime}_{1}(p,L)$, and so results of lattices smaller than $L=100$ were excluded in the fitting. This kind of deviation is also observed in BP of $m=3$ in the other kinds of lattices. The percolation threshold can also be found by the value of $p$ at the crossing points of $P_{w1}(p,L)$ and $P_{w2}(p,L)$ with various linear size ($L$) [28], as shown in Fig. 2. We ruled out the data from small lattices to reduce possible error from the finite-size-effect. All the values of the percolation threshold obtained by these methods were consistent with each other within error bars. The maximum of the probability of percolation only in one direction, $P_{w2}(p,L)-P_{w1}(p,L)$, which has negligible finite-size-effect, gives the most accurate percolation threshold.

The final estimate of the percolation threshold was obtained by taking an average. Table 1 shows the percolation threshold of all continuous transitions; cases with the first-order transition (the BP with $m>m_{c}$ and DP with $k<k_{c}=\Delta+1-m_{c}$) are omitted in the table. Note that $m_{c}=\Delta/2$ for even-coordinated lattices and $m_{c}=(\Delta+1)/2$ for odd-coordinated lattices, with exception of T7 (bounce lattice), which has $\Delta=4$ and $m_{c}=3$. The percolation threshold results of the BP with $m=1$ or $m=2$ are the same as that of the CP within error bars, as is expected. We confirmed that they are also consistent with exact or the most precise numerical results of the CP [23, 32, 41, 42] within relative errors of the order of $0.001\%$. The percolation threshold values of the BP and DP of other values of $m$ and $k$ for square, triangular, and honeycomb lattices are also consistent with references [6, 16, 17], but our work is much more precise.

Figure 4 shows maximum values of the derivative of percolation probability, $dP_{w1}(p,L)/dp$ and $dP_{w2}(p,L)/dp$, and derivative of percolation probability and strength of the largest cluster $P_{\infty}(p,L)$ at the percolation threshold $p_{c}^{(\infty)}$ calculated in this work. They satisfy the scaling relations [35, 29]:

[TABLE]

Therefore, the critical exponents can be obtained by fitting. In cases of the BP of $m=3$, there exists small but systematic deviation from the scaling behavior for $P_{\infty}(p_{c}^{(\infty)},L)$ in small lattices ($L<100$), which were excluded from the fitting. Table 1 summarizes the critical exponents obtained in this work. They are consistent with those of the two-dimensional classical percolation model ($\nu=4/3$ and $\beta/\nu=5/48$) within error bars, which are much smaller than references [10, 11, 13, 6, 15, 16, 17]. Therefore, we are convinced that continuous transitions of the BP and DP have the same critical exponents as the CP in two dimensions.

4 Summary

We extended the Newman-Ziff algorithm of the classical percolation to propose very efficient algorithms for the BP and DP models. Using these algorithms we studied the BP and DP in the eleven Archimedean lattices. The BP with $m\leq m_{c}$ and the DP with $k\geq(\Delta+1-m_{c})$ have continuous percolation phase transitions. We found that $m_{c}=\lfloor(\Delta+1)/2\rfloor$ except for the bounce lattice (T7), which has $m_{c}=\lfloor(\Delta+1)/2\rfloor+1$. Through the finite-size-scaling, we calculated the percolation threshold and critical exponents for the BP and DP with the continuous phase transition. We found that the critical exponents $\nu$ and $\beta$ are the same as those of the CP within error bars to conclude that the BP and DP almost certainly belong to the same universality class as the CP in two dimensions.

The algorithms presented in this paper can be directly applied to any dimensions and graphs. Since the BP and DP models are useful both in materials on lattices and in complex systems, studies of the BP and DP using these new algorithms in three-dimensional lattices and complex networks would be also interesting.

Acknowledgments

This work was supported by GIST Research Institute (GRI) grant funded by the GIST in 2019.