Percolation Threshold

Percolation theory studies connectivity in random networks. Applied to urban infrastructure, it answers a critical question: how many roads (or links) can fail before the city becomes disconnected? The percolation threshold marks the sharp transition between a connected and fragmented network.

1. Bond and Site Percolation

In bond percolation, each edge of a graph is independently open (functional) with probability $p$ or closed (failed) with probability $1 - p$. In site percolation, it is the nodes that are independently occupied with probability $p$.

The central quantity is the giant component -- the largest connected cluster of open edges/sites. Define $P_\infty(p)$ as the fraction of nodes in the giant component:

$$P_\infty(p) = \lim_{N \to \infty} \frac{|C_{\max}|}{N}$$

There exists a critical occupation probability $p_c$:

$$P_\infty(p) = \begin{cases} 0 & p < p_c \\ > 0 & p > p_c \end{cases}$$

For regular lattices, the critical thresholds are known exactly in some cases:$p_c^{\text{bond}} = 1/2$ for the square lattice (bond percolation),$p_c^{\text{site}} \approx 0.5927$ for the square lattice (site percolation), and $p_c^{\text{bond}} = 2\sin(\pi/18) \approx 0.3473$ for the triangular lattice.

2. Molloy-Reed Criterion

For random graphs with arbitrary degree distribution $P(k)$, the Molloy-Reed criterion determines when a giant component exists:

$$\frac{\langle k^2 \rangle}{\langle k \rangle} > 2 \qquad \Longleftrightarrow \qquad \text{giant component exists}$$

Derivation

Consider following a random edge to a node with degree $k$. The probability of reaching a degree-$k$ node is proportional to $k P(k)$ (biased sampling). The excess degree (number ofother edges) of this node is $k - 1$. The expected number of new nodes reachable from a randomly followed edge is:

$$\kappa = \frac{\sum_k k(k-1)P(k)}{\sum_k k P(k)} = \frac{\langle k^2 \rangle - \langle k \rangle}{\langle k \rangle}$$

For the giant component to grow, each step must reach at least one new node: $\kappa > 1$:

$$\frac{\langle k^2 \rangle - \langle k \rangle}{\langle k \rangle} > 1 \quad \Rightarrow \quad \frac{\langle k^2 \rangle}{\langle k \rangle} > 2$$

For Erdős-Rényi random graphs where each pair connects with probability $p$, we have $\langle k \rangle = (n-1)p$ and $\langle k^2 \rangle = \langle k \rangle + \langle k \rangle^2$. The criterion gives:

$$p_c = \frac{1}{\langle k \rangle} = \frac{1}{n - 1} \approx \frac{1}{n}$$

This is the classic result: in a random graph on $n$ nodes, the giant component appears when the average degree exceeds 1.

3. Critical Exponents and Scaling

Near the critical point, the percolation transition is characterised by power laws with universal critical exponents:

$$P_\infty(p) \sim (p - p_c)^{\beta} \qquad \text{for } p > p_c$$

$$\langle s \rangle \sim |p - p_c|^{-\gamma} \qquad \text{(mean cluster size)}$$

$$\xi \sim |p - p_c|^{-\nu} \qquad \text{(correlation length)}$$

In 2D (square lattice), the exponents are exactly known:

$\beta = 5/36$ (order parameter)
$\gamma = 43/18$ (susceptibility)
$\nu = 4/3$ (correlation length)

At the critical point itself, the cluster size distribution follows a power law:$\;n_s \sim s^{-\tau}$ with $\tau = 187/91 \approx 2.055$ in 2D. The fractal dimension of the percolating cluster at criticality is$d_f = 91/48 \approx 1.896$.

4. Urban Infrastructure Resilience

Percolation theory provides the framework for understanding urban infrastructure resilience. Key applications:

Road network failure: Bond percolation models random road closures (accidents, floods). The percolation threshold $p_c$ gives the fraction of roads that must remain open for the city to stay connected.
Targeted attacks: Removing high-betweenness nodes first is far more damaging than random failure. Scale-free networks are robust to random failure but vulnerable to targeted attack.
Cascading failures: When a road closes, traffic reroutes, potentially overloading adjacent roads -- triggering further failures. This cascade can push the system past $p_c$ even from a small initial perturbation.

For a street network with average degree $\bar{k} \approx 4$ (grid city), the random failure threshold is approximately $1 - p_c \approx 1 - 1/(\bar{k}-1) = 2/3$, meaning roughly two-thirds of roads can fail randomly before the network fragments. Suburban networks with $\bar{k} \approx 2.5$ are much more vulnerable: $1 - p_c \approx 1/3$.

5. Python: Bond Percolation Simulation

We simulate bond percolation on a square lattice, find the critical threshold, and plot the giant component fraction as a function of occupation probability.

Bond Percolation on Square Lattice

Python

script.py116 lines

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

np.random.seed(42)

# --- Bond percolation on square lattice ---
L = 40   # lattice size
G_full = nx.grid_2d_graph(L, L)
n_nodes = L * L
n_edges = G_full.number_of_edges()

p_values = np.linspace(0.0, 1.0, 50)
n_trials = 5

giant_fraction = np.zeros((len(p_values), n_trials))
second_fraction = np.zeros((len(p_values), n_trials))
n_clusters_arr = np.zeros((len(p_values), n_trials))

for trial in range(n_trials):
    for ip, p in enumerate(p_values):
        # Keep each edge with probability p
        G_perc = nx.Graph()
        G_perc.add_nodes_from(G_full.nodes())
        for e in G_full.edges():
            if np.random.random() < p:
                G_perc.add_edge(*e)

components = sorted(nx.connected_components(G_perc), key=len, reverse=True)
        giant_fraction[ip, trial] = len(components[0]) / n_nodes
        if len(components) > 1:
            second_fraction[ip, trial] = len(components[1]) / n_nodes
        n_clusters_arr[ip, trial] = len(components)

mean_giant = np.mean(giant_fraction, axis=1)
std_giant = np.std(giant_fraction, axis=1)
mean_second = np.mean(second_fraction, axis=1)
mean_clusters = np.mean(n_clusters_arr, axis=1)

# --- Find critical threshold (max of second-largest cluster) ---
pc_idx = np.argmax(mean_second)
pc_est = p_values[pc_idx]

# --- Plot ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Top left: Giant component fraction
ax = axes[0, 0]
ax.errorbar(p_values, mean_giant, yerr=std_giant, fmt='o-', color='#10b981',
            linewidth=2, markersize=4, capsize=3, markerfacecolor='#34d399')
ax.axvline(x=0.5, color='#f59e0b', linestyle='--', linewidth=1.5, label='$p_c = 0.5$ (exact)')
ax.axvline(x=pc_est, color='#ef4444', linestyle=':', linewidth=1.5, label=f'$p_c \\approx {pc_est:.2f}$ (est)')
ax.set_xlabel('Occupation probability p', fontsize=12, color='white')
ax.set_ylabel('$P_\\infty(p)$', fontsize=12, color='white')
ax.set_title('Giant Component Fraction', fontsize=13, color='white')
ax.legend(fontsize=10, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Top right: Second-largest cluster (peaks at p_c)
ax = axes[0, 1]
ax.plot(p_values, mean_second, 'o-', color='#06b6d4', linewidth=2, markersize=4)
ax.axvline(x=0.5, color='#f59e0b', linestyle='--', linewidth=1.5)
ax.set_xlabel('Occupation probability p', fontsize=12, color='white')
ax.set_ylabel('Second cluster / N', fontsize=12, color='white')
ax.set_title('Second-Largest Cluster (peaks at $p_c$)', fontsize=13, color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Bottom left: Number of clusters
ax = axes[1, 0]
ax.plot(p_values, mean_clusters, 'o-', color='#f59e0b', linewidth=2, markersize=4)
ax.axvline(x=0.5, color='#10b981', linestyle='--', linewidth=1.5)
ax.set_xlabel('Occupation probability p', fontsize=12, color='white')
ax.set_ylabel('Number of clusters', fontsize=12, color='white')
ax.set_title('Cluster Count', fontsize=13, color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Bottom right: Visualise percolation at p_c
ax = axes[1, 1]
G_crit = nx.Graph()
G_crit.add_nodes_from(G_full.nodes())
for e in G_full.edges():
    if np.random.random() < 0.5:
        G_crit.add_edge(*e)

components_crit = sorted(nx.connected_components(G_crit), key=len, reverse=True)
pos = {(i, j): (i, j) for i in range(L) for j in range(L)}

# Color largest component differently
node_colors = []
giant_set = components_crit[0]
for node in G_crit.nodes():
    if node in giant_set:
        node_colors.append('#10b981')
    else:
        node_colors.append('#334155')

nx.draw_networkx_edges(G_crit, pos, ax=ax, edge_color='#475569', width=0.3, alpha=0.4)
nx.draw_networkx_nodes(G_crit, pos, ax=ax, node_size=5, node_color=node_colors)
ax.set_title(f'Percolation at $p_c = 0.5$ (L={L})', fontsize=12, color='white')
ax.set_aspect('equal')
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)

fig.patch.set_facecolor('#0f172a')
for ax_row in axes:
    for ax in ax_row:
        ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()

print(f"Estimated p_c = {pc_est:.3f} (exact: 0.500)")
print(f"Giant component at p=0.6: {mean_giant[np.argmin(np.abs(p_values-0.6))]:.3f}")
print(f"Giant component at p=0.4: {mean_giant[np.argmin(np.abs(p_values-0.4))]:.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Large-Scale Percolation: Union-Find Algorithm

For large lattices, we need an efficient algorithm. The Union-Find (disjoint set) data structure processes percolation in near-linear time $O(N \alpha(N))$ where $\alpha$ is the inverse Ackermann function (effectively constant).

Large-Scale Percolation with Union-Find

Python

script.py118 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

class UnionFind:
    """Union-Find with path compression and union by rank."""
    def __init__(self, n):
        self.parent = list(range(n))
        self.rank = [0] * n
        self.size = [1] * n

def find(self, x):
        while self.parent[x] != x:
            self.parent[x] = self.parent[self.parent[x]]  # path compression
            x = self.parent[x]
        return x

def union(self, x, y):
        rx, ry = self.find(x), self.find(y)
        if rx == ry:
            return
        if self.rank[rx] < self.rank[ry]:
            rx, ry = ry, rx
        self.parent[ry] = rx
        self.size[rx] += self.size[ry]
        if self.rank[rx] == self.rank[ry]:
            self.rank[rx] += 1

def largest_component(self):
        return max(self.size[self.find(i)] for i in range(len(self.parent)))

def percolation_sweep(L, p_values, n_trials=3):
    """Efficient percolation using Union-Find."""
    N = L * L
    results = np.zeros((len(p_values), n_trials))

for trial in range(n_trials):
        for ip, p in enumerate(p_values):
            uf = UnionFind(N)

for i in range(L):
                for j in range(L):
                    node = i * L + j
                    # Right neighbour
                    if j < L - 1:
                        if np.random.random() < p:
                            uf.union(node, node + 1)
                    # Down neighbour
                    if i < L - 1:
                        if np.random.random() < p:
                            uf.union(node, node + L)

results[ip, trial] = uf.largest_component() / N

return np.mean(results, axis=1), np.std(results, axis=1)

# Sweep for multiple lattice sizes
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

p_values = np.linspace(0.3, 0.7, 40)
sizes = [20, 40, 80]
colors = ['#06b6d4', '#10b981', '#f59e0b']

ax = axes[0]
for L, color in zip(sizes, colors):
    mean, std = percolation_sweep(L, p_values, n_trials=3)
    ax.plot(p_values, mean, 'o-', color=color, linewidth=2, markersize=3, label=f'L = {L}')
    ax.fill_between(p_values, mean - std, mean + std, alpha=0.15, color=color)

ax.axvline(x=0.5, color='#ef4444', linestyle='--', linewidth=1.5, label='$p_c = 0.5$')
ax.set_xlabel('Bond probability p', fontsize=12, color='white')
ax.set_ylabel('Giant component fraction', fontsize=12, color='white')
ax.set_title('Finite-Size Scaling of Percolation', fontsize=13, color='white')
ax.legend(fontsize=10, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Power-law fit near p_c for largest L
ax2 = axes[1]
L_fit = 80
mean_fit, _ = percolation_sweep(L_fit, p_values, n_trials=3)

# Above p_c
mask = (p_values > 0.52) & (p_values < 0.68) & (mean_fit > 0.01)
p_above = p_values[mask]
P_above = mean_fit[mask]
log_dp = np.log(p_above - 0.5)
log_P = np.log(P_above)

# Linear fit in log-log
coeffs = np.polyfit(log_dp, log_P, 1)
beta_est = coeffs[0]

p_fit = np.linspace(0.51, 0.68, 50)
P_fit = np.exp(coeffs[1]) * (p_fit - 0.5)**beta_est

ax2.plot(p_above - 0.5, P_above, 'o', color='#10b981', markersize=6)
ax2.plot(p_fit - 0.5, P_fit, '--', color='#f59e0b', linewidth=2,
         label=f'$P_\\infty \\sim (p - p_c)^{{\\beta}}$, $\\beta \\approx {beta_est:.2f}$')
ax2.set_xlabel('$p - p_c$', fontsize=12, color='white')
ax2.set_ylabel('$P_\\infty$', fontsize=12, color='white')
ax2.set_title(f'Critical Scaling (L={L_fit})', fontsize=13, color='white')
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.legend(fontsize=11, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2, which='both')

fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()

print(f"Estimated critical exponent beta = {beta_est:.3f} (exact: {5/36:.3f})")
print(f"Union-Find complexity: O(N * alpha(N)) per trial")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Percolation on General Networks

For networks with degree distribution $P(k)$, bond percolation with probability $p$ effectively replaces $P(k)$ with a thinned distribution. The generating function approach gives the critical threshold:

$$p_c = \frac{\langle k \rangle}{\langle k^2 \rangle - \langle k \rangle} = \frac{1}{\kappa}$$

where $\kappa$ is the branching ratio from Section 2. For scale-free networks with $P(k) \sim k^{-\gamma}$ and $\gamma < 3$, the second moment diverges ($\langle k^2 \rangle \to \infty$), so:

$$p_c \to 0 \qquad \text{for scale-free networks with } \gamma < 3$$

This means scale-free networks are infinitely robust to random failure -- you can remove almost all edges and the giant component persists. However, targeted removal of the highest-degree nodes (hubs) quickly fragments the network. This has profound implications for urban infrastructure: transportation networks dominated by a few major highways are resilient to random disruptions but catastrophically vulnerable to targeted attacks on key junctions.

← Spectral Analysis Part VII →