Percolation on complex networks: Theory and application

Imagine a large porous stone immersed in water. Does the water come into the core of the stone? If so, there must be some paths formed by the pores running through the stone. However, as a whole, the connection of any two adjacent pores is probabilistic, supposing the probability is $p$. The problem thus reduces to whether there exist such paths for a given probability $p$. This actually is a typical example of percolation problems. The connection of two adjacent pores, in the terminology of percolation, is called occupying the corresponding bond between the two pores, and hence $p$ is the occupied probability. Obviously, if $p$ is large enough, the core of the porous stone can be wetted. In that case, the connected pores are able to form a cluster that penetrates the stone. Accordingly, percolation theory is mainly concerned with the existence of such a cluster and its structure with respect to the occupied probability $p$.

The above process is just a special case of the percolation process. The same question can also be proposed for many other systems constituted by random medium. Another typical example is the forest fire model. Suppose a burning tree can only ignite the trees in the adjacent sites, the destructiveness of the forest fire depends on the density of the trees, i.e., the probability $p$ of finding a tree on a site. This is also obviously a percolation process – when the adjacent trees form a cluster that can penetrate the entire forest, the fire could raze almost the entire forest; otherwise, the fire is constrained in a small area. Furthermore, similar regulations can also be applied to model the spreading of epidemics among individuals, where infected individuals can infect their neighbors probabilistically.

Although the soaking of porous stone, the forest fire and the spreading of epidemics seemingly belong to separated fields, the three are now converging in an intriguing manner. That is, is there a cluster of connected sites through the system? As is quite typical, it is actually easier to examine such an infinite cluster in an infinite system than just large ones. With the increasing of the probability $p$, there must be a critical $p_{c}$, called percolation threshold or critical point, below which such a cluster cannot be found. For a more detailed understanding of this criticality, the percolation model needs to be defined explicitly. That is what we are going to talk about.

To facilitate presentation, we consider a two-dimensional lattice here as shown in Fig.1 (a). Mathematically, each bond is occupied with probability $p$ or unoccupied with probability $1-p$. Then, the occupied bonds connect the sites into clusters. This model is called bond percolation, which can be used to model the process of the soaking of porous stone and the spreading of epidemics. For the forest fire model, one often uses a slightly different percolation model, the so-called site percolation model. For this model, we occupy each site with probability $p$ rather than bonds as shown in Fig.1 (b). In general, the bond percolation is considered less general than the site percolation due to the fact that the bond percolation can be reformulated as a site percolation on a different lattice, but not vice versa.

The percolation theory mainly focuses on the emergence of the infinite cluster with the increasing of probability $p$. To characterize this phenomenon, one often adopts the size of the giant cluster, which is defined as

Here, $N$ is the size of the system (site number), and $s_{1}$ is the number of sites in the largest cluster. As shown in Fig.1, with the increasing of $p$, there must be a critical point $p_{c}$, above which a non-zero $S$ can be found. This figures out the percolation transition of the system with respect to the control parameter $p$, and $S$ is the corresponding order parameter.

In addition, there is another commonly used order parameter called wrapping probability, which is defined as the probability that a cluster wraps around the periodic boundary conditions on a regular lattice. For large systems, this probability is equal to the probability that the system percolates. This parameter is usually used to estimate the position of the percolation threshold, since it is a size-independent parameter at the critical point. Specifically, cluster wrapping can be defined in a number of different ways, such as wrapping around one direction, wrapping around one direction but not the other, and wrapping around both directions. This probability, however, cannot be well defined in network systems without spatial constraints, so that the size of the giant cluster is the preferred order parameter in network science.

To describe the features of finite clusters, the distribution of cluster sizes $p_{s}=m_{s}/\sum_{s}m_{s}$ is also used, where $m_{s}$ is the number of the clusters with size $s$. Sometimes, the normalized cluster number $n_{s}=m_{s}/N$ is also used to feature the cluster size distribution. It is obvious that $p_{s}=n_{s}N/\sum_{s}m_{s}=n_{s}\langle s\rangle$, where $\langle s\rangle=N/\sum_{s}m_{s}=\sum_{s}sp_{s}$ is the average cluster size.

It must be pointed out that the average size $\langle s\rangle$ is not the mean cluster size $\chi$ commonly used in percolation theory, which is defined as $\chi=\sum_{s}s^{2}n_{s}$. It is not hard to know that $sn_{s}$ is the probability that a randomly chosen site belongs to a cluster with size $s$. Thus, the mean cluster size $\chi$ actually is the average size of the cluster that a randomly chosen site belongs to. Together with the characteristic length $\xi$ and the characteristic size $s_{\xi}$, above which the clusters are exponentially scarce, these statistics are often used to describe and characterize the percolation transition [7]. The critical behaviors of these parameters will be briefly introduced in the next section.

Besides, there exist several other variants of percolation. For example, one could drop the assumption of independence of occupation, so that the occupation of a site or a bond could depend on the state of other sites or bonds. This type of percolation is often called dependent percolation, which is widely used in network science to model the spreading dynamics and the cascading process [3, 5, 6, 8, 23]. Such models will also be discussed in subsequent parts of this paper.

The prerequisite for the study of critical phenomena is determining the percolation threshold $p_{c}$. Figure 1 indicates that the site percolation and the bond percolation have different critical points. In fact, the threshold of bond percolation on square lattice can be solved exactly by the duality transformation or real-space renormalization [24, 7, 25, 26, 27], that is $p_{c}=1/2$. However, there is no exact solution for site percolation on square lattice. Through Monte Carlo simulation, a good estimate of the critical point can be obtained, such as $p_{c}=0.59274621(13)$ in Ref.[28], and $p_{c}=0.5927465(4)$ in Ref.[29]. In Tab.1, we list the percolation thresholds for some typical systems.

In the subcritical regime ($p<p_{c}$), all clusters are finite and the size distribution has a tail which decays exponentially; in the supercritical regime ($p>p_{c}$), an infinite cluster can be found in the system, and the size distribution of other finite clusters has a tail which decays slower than exponential. Near the critical point, some asymptotic behaviors can be found, referred to as critical phenomena. In the percolation theory, these behaviors are characterized by the critical exponents [7],

with relations

Here, $\beta$, $\gamma$, $\nu$, $\sigma$, and $\tau$ are the so-called critical exponents, which determine the universal class of the percolation transition. Besides, at the critical point the characteristic length $\xi$ and size $s_{\xi}$ also have a relation

The exponent $d_{f}$ is often called fractal dimension [7, 49], characterizing the structure of the infinite cluster at the critical point. Assuming the dimension of the system is $d$, there is another relation between critical exponents, called hyperscaling,

Note that this relation is also universal, and independent of the topological structure of the system. The phase transition theory points out that there is an upper critical dimension $d_{c}$ (for percolation $d_{c}=6$), above which the critical exponents become the same as that in mean-field theory [50, 51]. So this hyperscaling relation holds only for $d\leq d_{c}$. It is also worth mentioning that the geometric structure of high-dimensional percolation clusters cannot be fully accounted for by the counterparts of random networks, though both of them have the mean-field nature [52].

In addition, the critical behavior below $d_{c}$ is different from the mean-field approximation which is valid away from the phase transition. For these systems, the renormalization group theory has made remarkable predictions about the behavior of the percolation near and at the threshold [53, 54]. In the renormalization group theory, there is also a lower critical dimension (for percolation $d_{l}=2$), below which there is no phase transition.

Furthermore, the universality principle states that the value of $p_{c}$ is determined by the local structure of the system, whereas the behavior of clusters that is observed near $p_{c}$ is independent of the local structure (lattice type and percolation type). In this sense, percolation is believed to be a substrate-dependent but model-independent process, and therefore, the critical exponents of the transition are only determined by the geometry of the system, and identical for the bond and site percolation. The critical exponents are thus more natural to be considered than the threshold $p_{c}$, and there is no need to deal with the site and bond percolation, individually.

To end this subsection, it must be pointed out that the critical exponent $\beta$ is distinguishable for the site and bond percolation on scaling-free (SF) networks [55]. This is a special case derived from the vanished percolation threshold, which has no effect on the other properties discussed above.

In network science, the elements of a system and the connections between them are no longer known as site and bond. Instead, they are often called node and link, or vertex and edge, respectively. The number of links a node has is called degree, labeled $k$.

To exactly represent the connection pattern of a network, one often uses a $N\times N$ matrix $\mathbf{A}$, called adjacent matrix, where $N$ is the number of nodes in the network. In this way, the element $a_{ij}$ of $\mathbf{A}$ is one when there is a link from node $i$ to node $j$, and zero when there is no link. For undirected networks, it must have $a_{ij}=a_{ji}$, thus $\mathbf{A}$ is a symmetric matrix. For weighted networks, the element $a_{ij}$ can further be any non-zero value to represent the corresponding link weight.

Instead of studying a special network with a given adjacent matrix, the network science is more on discovering the common nature of a class of networks, i.e., a network ensemble. For that, networks are often featured by the degree distribution $p_{k}$, which provides the probability that a randomly selected node in the network has degree $k$.

A typical network ensemble is the one with Poisson degree distribution

where $\langle\rangle$ means the average over all the nodes. This just is the ensemble of the known Erdős–Rényi (ER) random networks/graphs [47], and can be simply realized by randomly connecting a given number of links among a set of nodes, or connecting each pairs of nodes with a given probability. Quite obviously, the generation of ER networks is actually a percolation process. Only when the system percolates ($\langle k\rangle>1$), the giant cluster can be found. Specifically, when $\langle k\rangle<1$ almost surely each cluster has size $\log N$, at $\langle k\rangle=1$, the largest cluster has a size of order $N^{2/3}$, and at $\langle k\rangle>1$, there exists a single giant cluster of size of order $N$ and all other clusters have a size of order $\log N$. Since such systems have no spatial constraint, the critical phenomena in this ensemble just recover the mean-field solution.

Besides, distinguishing from the degree distribution, there is another kind of commonly used networks, which represents a deeper organizing principle of real networks called the SF property [2]. In mathematical terms, the SF property translates into a power-law function of the form

For a real network, there must also be an upper bound of degree, called degree cutoff $K$.

The main difference between an ER and a SF network comes in the tail of the degree distribution. For ER networks, most nodes have comparable degrees and hence the degree cutoff is in the order of the average degree. In contrast, nodes with very large degrees are expected in SF networks, called hubs of the networks. Indeed, the degrees of hubs grow with the network size, and thus can grow quite large. This indicates that SF networks have strong heterogeneity, and hence the percolation transition bears a strong dependence upon the degree distribution. Beyond this, there are many other measurements to further subdivide the network ensembles, such as clustering, correlation, and community. In recent years, temporal networks, multiplex networks, and high-order networks have also received a lot of attention. Because of space limitation and the focus of this paper, the details will not be included here and brief discussions will be given in this article when necessary. Further details on these can also be found in Refs.[3, 5, 6, 1, 2].

After a review of the fundamental properties of complex networks, we explore the mathematical modeling of network ensembles, and mainly focus on the configuration model and the hidden parameter model [1, 2]. In principle, the two models can generate any network ensembles with a meaningful degree distribution. In fact, these two models are usually only used to realize SF property, since there exist easy-implemented models for Poisson degree distribution (ER networks) and small-world (SW) networks [1, 2]. Note that the SF network also has an easy-implemented model, Barabási–Albert (BA) model [56], nevertheless, the generated SF network has its own structural features and limitations rather than a network with programmable and tunable SF properties.

From random graph theory [47], there are two ways to generate a network with Poisson degree distribution. One is randomly connecting nodes until the excepted number of links is met, and the other is connecting each pair of nodes with a given probability. For a given average degree, networks implemented by the former method must be a subset of those by the later one. In spite of this, the two methods will become equivalent in the thermodynamic limitation ($N\to\infty$).

The SW network grows out of a regular network, such as the triangle lattice and the square lattice, by randomly rewiring a small fraction of links [57, 58, 59, 60]. Since there is no spatial constraint for the rewired link, the average distance of nodes becomes much smaller than that of the underlying regular network. This is why it is called small-world. In turn, the regular backbone of the SW network usually leads to a high clustering, namely, the neighbors of a node are also connected. This is another characteristic of the SW network. For this type of networks, the degree distribution is often not the concern, and one often uses the backbone network and the rewiring probability to define a SW network ensemble.

The configuration model can help us build a network with a pre-defined degree distribution [61, 62, 1, 2]. The algorithm consists of two steps. First, assigning degrees to each node drawn randomly from the pre-defined degree distribution, represented as stubs. Then, two stubs are selected randomly and connected. Repeating this procedure until all stubs are paired up. If there is nothing in this procedure to forbid self- and multi-connections, the obtained network is probably not a simple graph as usually studied in network science. While rejecting such connections could add much more overhead to the program’s execution time.

Indeed, one of the effective methods to generate random networks by the configuration model without self- and multi-connections can be as follows. First, connecting the stubs in any fast and easy-to-implement ways (the details are dependent on the specific degree distribution), and recording the self- and multi-connections. Then, rewiring a fraction of links of the obtained network, and the self- and multi-connections must be included. Note that the rewiring leading to self- or multi-connections is forbidden.

A schematic diagram of the rewiring procedure is shown in Fig.2. For the degree distribution with the maximum degree $K$ larger than $\sqrt{N}$, this method could be much more effective than the original one, which may not even complete the network. This is because there must be some degree correlation in the network with $K>\sqrt{N}$ [63, 64], randomly connecting stubs with self- and multi-connection forbidden is a very inefficient procedure. Of course, the efficiency of the algorithm also depends on the implementing way and the detailed property of the network, thus this just is a general discussion and not always the case.

By the way, the rewiring operation as shown in Fig.2 is often called degree preserving randomization, which can randomize the connections of a network without changing the degree distribution. With this operation, one can eliminate other properties from the network, for example, clustering, and assure that a certain network feature is predicted by its degree distribution alone. If the network is large enough and all the degrees are much smaller than the network size $N$, the degree preserving randomization could turn the network to be locally tree-like.

Note that the correlation derived from the degree distribution cannot be undone by this randomization procedure. A typical example is a network with $K>\sqrt{N}$. Consequently, when we study networks with hubs, i.e., SF networks, the configuration model and the degree preserving randomization must be carefully dealt with. In turn, with some constraints, this rewiring operation can also “randomize” the network into the one with a given high-order property [65, 66].

There is another model, hidden parameter model, to generate networks with a pre-defined degree distribution [67, 68, 69, 70, 71]. In this model, each node $i$ is assigned a hidden parameter $\eta_{i}$, chosen from a pre-defined distribution $p_{\eta}$. Then, we connect each node pair with probability $p(\eta_{i},\eta_{j})=\eta_{i}\eta_{j}/N\langle\eta\rangle$. The expected number of links is thus $L=\sum_{i,j}\eta_{i}\eta_{j}/2N\langle\eta\rangle$. For SF networks, we can set $\eta_{i}\propto i^{-\alpha}$, and the obtained network has the degree distribution $p_{k}\propto k^{-(1+1/\alpha)}$.

The SF networks generated by the configuration model and the hidden parameter model also have their own properties, thus are not exactly equivalent, at least for small systems. First of all, the boundaries of degrees must be given for the configuration model to normalize the degree distribution. However, the boundaries of degrees are not needed for the hidden parameter model, in which a degree cutoff occurs naturally at $N^{1/(\lambda-1)}$, thus known as the natural cutoff of the SF network [72, 73, 74, 75, 76]. This cutoff can be also found in the configuration model, if the upper boundary $M$ used to normalize the degree distribution is larger than $N^{1/(\lambda-1)}$. Then, the hidden parameter model generates no self- and multi-connections, since only a single connection would occur between each node pair. Another is that the degree distribution of the SF network generated by the hidden parameter model has negative deviations for small degrees, while the one generated by the configuration model almost perfectly matches the pre-defined degree distribution. In addition, setting the total number of links is not allowed in the configuration model, since it is fixed by the pre-defined degree distribution. Rather, the hidden parameter model allows controlling the link number, accurately. Finally, what need to be pointed out is that both the networks generated by the configuration model and by the hidden parameter model are only a possible realization of the pre-defined degree distribution. To check the property induced by a certain degree distribution, the ensemble average must be done.

In recent years many works have also been devoted to the so-called multiplex networks (or multilayer networks, interdependent networks) [8, 9]. In this network ensemble, links are divided into different layers to represent different types of connections between nodes. To generate such a network ensemble, one can first generate several networks using the methods mentioned above, and then bundle nodes from different networks together to form the layered connections. It can also be interpreted as inserting some inter-connections between nodes in different networks (layers), which represent the dependence relations between them. The bundling of nodes (or inter-connection) could be one-to-one, one-to-many, or with some other rules, that depending on the inter-correlation between layers. Of course, there are many other network models with different properties, one can find them in special books and reviews about complex networks.

The percolation process is tantamount to dilute links and thus changes the degree distribution $p_{k}$. So, for convenience and generality, the occupation of nodes or links will not be considered in the following discussion. After obtaining the result, we can revise it by considering the diluting effect of the percolation on the degree distribution $p_{k}$.

When a network percolates, there must be an infinite cluster, in which the branching process is endless. Specifically, when arriving at a node by following a link, the node must have some other links (at least one) to ensure the continuity of branching. In terms of this, assuming a link belongs (connects) to the giant cluster with probability $R$, we have a self-consistent equation

where $q_{k}=p_{k}k/\langle k\rangle$ is the probability that the node reached by following a link has degree $k$, known as the excess-degree distribution, and $G_{1}(x)=\sum_{k}q_{k}x^{k-1}$ is the corresponding generating function. The term $1-(1-R)^{k-1}$ in the square brackets just means that at least one excess link is needed to keep branching.

After obtaining $R$ from Eq.(13), the order parameter $S$, i.e., the probability that a randomly chosen node belongs to the giant cluster (the infinite cluster), can be expressed as

where $G_{0}(x)=\sum_{k}p_{k}x^{k}$ is the generating function of the degree distribution. The right hand side of this equation means that at least one of the neighbors must belong to the giant cluster. From the perspective of branching, this equation can also be interpreted as the start point of the branching process, while Eq.(13) is the intermediate process. In the infinite cluster any nodes can be the start point, so $S$ is an average probability, as well as the parameter $R$.

Furthermore, for the percolation with occupied probability $p$, we only need to replace the generating functions in Eqs.(13) and (14) with those of the diluted network. Next, let us calculate the degree distribution of the diluted network. Regardless of the type of the percolation (site percolation or bond percolation), the diluting ratio of links is $p$, i.e., each link is preserved with probability $p$. So the generating functions $g_{0}(x)$ and $g_{1}(x)$ of the diluted networks can be written as

Then, inserting them into Eqs.(13) and (14), we have [48, 87]

Before continuing to consider the details of these two equations, we must point out that the order parameter $S$ here is the size respected to the diluted network. For bond percolation, the dilution is only for links, so the diluted network has the same size as the original network and $S$ is also the size respecting to the original network. However, for site percolation the diluted network consists of only occupied nodes (fraction $p$ respected to the original network). As a result, Eq.(18) should be revised for site percolation, that is

In general, one can solve Eq.(17) first, and then insert $R$ into Eq.(18) or (19) to get the order parameter $S$ for the bond percolation and the site percolation, respectively. From Eqs.(17)-(19), we can also find that the two types of percolation models give the same probability $R$, but different giant clusters. For mathematical tractability, sometimes one just uses $R$ to represent the meaning of $pR$ for site percolation in Eqs.(17) and (19), that is

These two equations are fully equivalent to Eqs.(17) and (19).

Although the mean-field approximation is used in the above discussion, i.e., all the nodes have the same probability $S$ and all the links have the same probability $R$, Eqs.(17)-(21) can give an exact solution for the percolation on tree-like networks, see Fig.3 as an example. A comparison between the mean-field prediction and the numerical simulation on some real-world networks can also be found in Ref.[88].

Following the idea of branching, the microstructure of the giant cluster in random networks [89], such as degree distribution, and degree correlations, as well as that of temporal networks [90], can also be obtained. For networks with loops, different links could lead to the same nodes in the branching process. In other words, the branching processes starting from different links are no longer independent of each other, thus all the equations (13)-(21) are not valid. That is why the above method is only applicable to the tree-like structure.

Besides the giant cluster, the percolation theory also concerns the behaviors of finite clusters, such as mean cluster size, and cluster size distribution. Next, along with the idea of the branching process used above, we will further show how to obtain the information of finite clusters in percolation model.

As discussed in Ref.[48], it is more convenient to study the distribution $\pi_{s}=n_{s}s=p_{s}s/\langle s\rangle$ rather than $n_{s}$ or $p_{s}$, which means the probability distribution of the size of the cluster that a randomly chosen node belongs to. One can further define $\rho_{s}$ as the size distribution of the cluster at the end of a link. Accordingly, the generating functions have the forms $H_{0}(x)=\sum_{s=0}\pi_{s}x^{s}$ and $H_{1}(x)=\sum_{s=0}\rho_{s}x^{s}$. For convenience, the occupation of links and nodes is not taken into account in the following. When necessary, one can use the revised generating functions Eqs.(15) and (16) to get the corresponding results.

Note that a zero size has no physical senses, so $\pi_{0}=0$, while $\rho_{0}$ can be non-zero corresponding a node with degree one. Since $H_{0}(1)$ and $H_{1}(x)$ only contain the finite clusters, and the giant cluster is excluded (if there is one), we have $H_{0}(1)=\sum_{s}\pi_{s}=1-S$ and $H_{1}(1)=\sum_{s}\rho_{s}=1-R$, where $S$ and $R$ are defined as that for Eqs.(13) and (14). Below the critical point, $R=0$ and $S=0$, so $H_{0}(1)=1$ and $H_{1}(1)=1$.

As the schematic in Fig.4, the cluster size distributions $\pi_{s}$ and $\rho_{s}$ are dependent on both the degree/excess-degree distribution and the branching size at the end of a link. Thus, Fig.4 can be translated as the equations

Note that there is one more factor $x$ on the right hand sides of the two equations for the contribution of the root node. Using the relations $H_{0}(1)=1-S$ and $H_{1}(1)=1-R$, these two equations just recover Eqs.(13) and (14) when $x=1$.

In principle, we can solve Eqs.(22) and (23) together to find $H_{0}(x)$, and then use

to find the cluster size distribution. However, usually no closed form can be found for $H_{0}(x)$. Instead, we can substitute Eqs.(22) and (23) into Eq.(24), then do a change of variables via Cauchy formula, and finally get

With this, the cluster size can be obtained directly from $G_{1}(x)$, therefore avoiding solve Eqs.(22) and (23) [91, 92, 93]. As $\pi_{s}\propto sn_{s}\propto sp_{s}$, the result must have the form $\pi_{s}\propto s^{1-\tau}$, where $\tau$ is the Fisher exponent.

Furthermore, the mean cluster size $\chi$ can also be obtained from Eqs.(22) and (23). From the definition, the mean cluster size $\chi$ can be represented as

In the last equation the differential of Eq.(22) is used. We can further employ Eq.(23) to replace the term $H_{1}^{\prime}(1)$, and yield

It is obvious that $\chi$ is divergent when

This is the Molloy–Reed criterion $\langle k(k-2)\rangle=0$ found in random graph theory [77], above which ($G_{1}^{\prime}(1)>1$) the giant cluster exists.

In addition, the discussion of finite clusters based on generating functions $H_{0}(x)$ and $H_{1}(x)$ can also work without the help of generating functions $G_{0}(x)$ and $G_{1}(x)$, so it is not limited to uncorrelated tree-like networks, see Sec.3.4 for an example of the percolation on growing networks.

It is known that the $q\to 1$ limit of the $q$-state Potts model without external field corresponds to the percolation problem, which is thus applied to the study of percolation on networks [95, 96, 16, 24]. To consider the percolation on a network $\mathcal{G}$, we introduce a $q$-state Potts model with the Hamiltonian

where $\alpha_{i}=1,2,3,\ldots,q$ is the spin state of node $i$, and the sums $\sum_{\langle i,j\rangle}$ and $\sum_{i}$ are over all the links and nodes of the network, respectively. $\delta_{x,y}$ is the Kronecker delta function, i.e., $\delta_{x,y}=0,1$ if $x\neq y$ and $x=y$, respectively. Here, the ferromagnetic interaction is used $J>0$, and the magnetic field $H>0$ is applied to the spin $\alpha_{H}$. Then, the partition function can be written as

where $K=\beta J$, $L=\beta H$, and the sum $\sum_{\{\alpha\}}$ is over all the spin configurations of the system. We can further expand the products and use the subnetworks of $\mathcal{G}$ to represent the terms in the expansion, which is often called Fortuin–Kasteleyn cluster representation [97].

Next, we give a brief derivation for the Fortuin–Kasteleyn cluster representation. For arguments $a_{i}$, $i\in I$, the product $\prod_{i\in I}(1+a_{i})$ is equivalent to the sum of all the possible polynomials formed by $a_{i}$, $i\in I$, so that $\prod_{i\in I}(1+a_{i})=\sum_{I^{\prime}\subseteq I}\prod_{i\in I^{\prime}}a_{i}$, where $I^{\prime}$ is a subset of $I$, and the sum $\sum_{I^{\prime}}$ is over all the possible cases. By this formula, the product $\prod_{\langle i,j\rangle}$ in Eq.(30) can be represented by the sum of the subnetworks $\mathcal{G}^{\prime}$ of $\mathcal{G}$, that is

where $l(\mathcal{G}^{\prime})$ is the number of links in subnetwork $\mathcal{G}^{\prime}$. Note that the term $\prod_{(i,j)\in\mathcal{G}^{\prime}}\delta_{\alpha_{i},\alpha_{j}}$ gives nonzero value only when the nodes in the same clusters are in the same states (any of the $q$ states $\alpha=1,2,\ldots,q$). This can further simplify the partition function to be

Here, the product $\prod_{c}$ is over all the clusters formed by the links in subnetwork $\mathcal{G}^{\prime}$, and $s_{c}$ is the number of nodes in cluster $c$. The sum over spin $\{\alpha\}$ in Eq.(30) has now been replaced by a sum over subnetwork $\{\mathcal{G}^{\prime}\}$, this is just the Fortuin–Kasteleyn cluster representation of the partition function of Potts model.

Let $e^{K}-1=p/(1-p)$, the partition function Eq.(32) can be rescaled as

where $l(\mathcal{G})$ is the total number of links in $\mathcal{G}$. If $p$ is interpreted as the occupied probability of links, the term $P_{\mathcal{G}^{\prime}}=p^{l(\mathcal{G}^{\prime})}(1-p)^{l(\mathcal{G})-l(\mathcal{G}^{\prime})}$ is just the probability of a percolation configuration with $l(\mathcal{G}^{\prime})$ occupied links. Thus, Eq.(33) bridges the Potts model and the percolation model.

From Eq.(33), we can write the free energy per node as

For convenience, we further define

For $L>0$ the sum in this equation is over all the finite clusters, then the size of the percolating cluster is

and the mean cluster size is

By the above two equations, the percolation property can then be extracted from the Potts model.

The product $\prod_{c}$ in Eq.(33) is over all the clusters. It thus can be rewritten in another form of

where the product $\prod_{s}$ is over the size of the clusters, and $m_{s}$ is the number of clusters with size $s$. By this formula, function $h(K,L)$ can also be expressed as

Here, $n_{s}=\langle m_{s}/N\rangle$ is the cluster size distribution. This is just the generating function $h(K,x)=\sum_{s}n_{s}x^{s}$ with $x=e^{-L}$, which can be used to study the cluster size distribution.

Specifically, once the Hamiltonian (Eq.(29)) of a network is known, all the percolation parameters can be obtained by the above equations. However, for most network models, it is difficult to simplify the sum over links in Eq.(29). A solvable case is the hidden parameter model, for which the pre-defined connection probability is given for each pair of nodes. Thus, the sum over links can be translated as a weighed (connection probability) sum over all the node pairs. With this technique, the Potts model formulation has been used in solving the percolation on SF networks [98, 99], on correlated hypergraphs [100], and in generalized canonical random network ensembles [101].

Through the so-called recurrence relation, the mean-field equations for tree-like networks can be recovered by the Potts model formulation [102]. For this purpose, the partition function can be represented as an integration of a root, labeled node $0$, and the branchings from it, that is

Here, $z_{i}(\alpha_{0})$ is the partition function of the branching from node $i$, and the product $\prod_{i}$ runs over all the root’s nearest neighbors (nodes in the first shell of node $0$). Due to the symmetry, $\alpha$ in the second term of Eq.(40) can be any state excluding $a_{H}$. Obviously, this formulation is only for the tree-like networks, if not, the product $\prod_{i}$ must include some double counting.

Different from Eq.(30), the partition function equation (40) represents the effective state of the root node. In this way, the free energy per nodes is $\ln Z$, then we have

To study the percolation model, we assume here that $z_{i}(\alpha)$ describes the state for $q=1$ and $L=0$, thus it is irrelevant to the differentials $\partial/\partial q$ and $\partial/\partial L$. For convenience, we further rescale $z_{i}(\alpha)$ by $z_{i}(\alpha_{H})$, i.e., $x_{i}\equiv z_{i}(\alpha)/z_{i}(\alpha_{H})$. Then, by Eq.(36), we can find the order parameter of the percolation transition

From this equation, the physical meaning of $x_{i}$ becomes clear, that is the probability that the branching from node $i$ is not the giant cluster.

To obtain the order parameter, we need to further find $x_{i}(\alpha)$. Comparing Eqs.(40) and (30), we have

where the sums $\sum_{j}$ and $\sum_{\langle j,k\rangle}$ run over the nodes and the links in the branching from node $i$, respectively. Note that the interaction between node $i$ and the root $K\delta_{\alpha_{0},\alpha_{i}}$ is also included in this partition function. Furthermore, it is not hard to find that $z_{i}(\alpha_{0})$ satisfies the following recurrence relation