---
title: Invaded-Cluster Algorithm in Statistical Physics
url: https://www.emergentmind.com/topics/invaded-cluster-algorithm
type: topic
---

# Invaded-Cluster Algorithm in Statistical Physics

The invaded-cluster algorithm is a Monte Carlo cluster method introduced by Machta–Chayes–Redner–Chayes for the Potts and Ising models on a torus and subsequently extended in two distinct directions: an equilibriumlike formulation for extracting critical properties in the canonical ensemble, and a homological generalization for $\mathbb Z_q$ Potts lattice gauge theory on the $4$-torus. Its defining feature is the replacement of a fixed bond-occupation step by a stopping rule tied to percolation: in the classical setting the algorithm grows a random subgraph until a nontrivial winding appears, whereas in the lattice-gauge setting it grows a plaquette subcomplex until a full basis of nontrivial $2$-cycles appears. In both forms, the method self-drives to criticality and performs nonlocal updates that strongly reduce critical slowing-down relative to single-spin or single-edge dynamics [1004.1509] [2507.13503].

## 1. Classical formulation in the Fortuin–Kasteleyn representation

For the $q$-state Potts model, the starting point is the Hamiltonian
$$
H=-J\sum_{\langle i,j\rangle}\bigl(\delta_{\sigma_i,\sigma_j}-1\bigr),\quad \sigma_i=1,\dots,q,\quad p=1-e^{-\beta J}.
$$
Under the Fortuin–Kasteleyn (FK) mapping, the partition function is rewritten as
$$
Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},
$$
where $B$ is the total number of edges, $b(\gamma)$ is the number of occupied bonds, and $c(\gamma)$ is the number of connected components. This is the representation that underlies standard cluster updates such as Swendsen–Wang (SW) [1004.1509].

In SW dynamics one alternates a bond-placement step and a cluster-flip step. For each nearest-neighbor pair $(i,j)$ with $\sigma_i=\sigma_j$, a bond is placed with probability $p$; then each FK-cluster is assigned an independent uniform spin label. The invaded-cluster (IC) algorithm modifies only the first part of this procedure. Instead of using a fixed $p$, it erases all bonds, identifies the satisfied edges,
$$
n_{ss}=\#\{\langle i,j\rangle:\sigma_i=\sigma_j\},
$$
and then places bonds one by one on uniformly random satisfied edges not yet bonded until a cluster spans the system, for example by wrapping around under periodic boundary conditions. The instantaneous estimator
$$
p_{\{\sigma\}}=\frac{b}{n_{ss}}
$$
is then recorded, the bonds are erased, SW-clusters are rebuilt, and the clusters are flipped at random [1004.1509].

This formulation makes the IC algorithm a self-adjusting variant of FK cluster dynamics. The stopping criterion is not prescribed by an external control parameter but by the onset of percolation in the current spin background.

## 2. Self-driving criticality and the non-equilibrium ensemble

The central motivation for the classical IC construction is that it drives the system to criticality without prior knowledge of the critical bond probability. After many steps, the average $\overline{p}_{\{\sigma\}}$ converges to $p_c(L)$, the finite-size critical bond probability. Because each update is forced to stop at percolation, the method “self-drives” to the critical point [1004.1509].

The standard IC algorithm also inherits the main advantage of cluster methods over local dynamics: the update is nonlocal. Rather than changing one degree of freedom at a time, it reconstructs a large correlated structure and then resamples it through the cluster-flip stage. This is the basis for its strong suppression of critical slowing-down in the original Potts-model setting and in later generalizations.

A central limitation, however, is that the induced ensemble is not canonical. The distribution of the instantaneous $p_{\{\sigma\}}$ is far wider than the canonical binomial variance $\sim 1/\sqrt{n_{ss}}$. As a consequence, the resulting ensemble is non-equilibrium and yields wrong finite-size scaling exponents [1004.1509]. This point is the main methodological caveat associated with the original IC algorithm: self-tuning to criticality does not by itself imply equilibrium sampling.

## 3. Equilibriumlike invaded-cluster algorithm

The equilibriumlike invaded cluster (EIC) algorithm was introduced precisely to retain self-driving while restoring the equilibrium ensemble. The construction groups Monte Carlo steps into blocks of $N_a$ successive IC-style moves and constrains the fluctuations of $p_{\{\sigma\}}$ to a window of size $2v$ around the previous block’s mean $\bar p_{i-1}$, with
$$
v=\tilde v\,L^{-d/2}.
$$
At the end of block $i$ one computes
$$
\bar p_i=\frac{1}{N_a}\sum_{\rm block}p_{\{\sigma\}},
$$
which becomes the center of the window for the next block [1004.1509].

The rationale is explicit. By choosing $v\sim L^{-d/2}$, one enforces
$$
\mathrm{var}(p_{\{\sigma\}})\approx L^{-d},
$$
matching the equilibrium binomial noise $\propto 1/n_{ss}\propto L^{-d}$. The self-driving mechanism ensures $\bar p_i\to p_c(L)$, while the fluctuation constraint ensures that the sampled ensemble is canonical. In this sense, EIC is not a different notion of criticality detection; it is a restriction on the allowed bond-fraction fluctuations so that percolation-based self-tuning is compatible with equilibrium finite-size scaling [1004.1509].

The critical-scaling results reported for EIC reflect this restoration of equilibrium behavior. For the largest-cluster mass,
$$
\overline{s_{\max}(L)}=a_s\,L^{y_h}+b_s\,L^{y_h+\omega_1}+\cdots,
$$
the fits give $y_h=2.4815(5)$ for the $3$D Ising case with $L\le 96$, and $y_h=1.8661(7)$ for the $2$D $q=3$ Potts case with $L\le 1024$. For the thermal exponent obtained from
$$
\frac{\partial\ln m}{\partial\beta}\sim a_{em}\,L^{1/\nu}+b_{em}\,L^{1/\nu+\omega_4},
$$
the fits give $1/\nu=1.586(5)$ in $3$D Ising and $1/\nu=1.201(8)$ in $2$D $q=3$ Potts. The corresponding three-parameter fits for the finite-size shift,
$$
p_c(L)=p_c(\infty)+a_p\,L^{-1/\nu}+\cdots,
$$
yield $p_c(\infty)=0.358097(1)$ for $3$D Ising and $p_c(\infty)=0.633975(1)$ for $2$D $q=3$. The reported exponents agree with known exact or best-numerical values to within $\sim 0.1\%$-$0.5\%$ [1004.1509].

## 4. Dynamical behavior, autocorrelations, and parameter dependence

The dynamical analysis of EIC is formulated in terms of the normalized autocorrelation function
$$
\Gamma_O(t)=\frac{\langle O(0)O(t)\rangle-\langle O\rangle^2}{\langle O^2\rangle-\langle O\rangle^2}.
$$
For $t\ll N_a$, EIC exhibits a single-exponential decay,
$$
\Gamma_O(t)\sim e^{-t/\tau_O},
$$
while around $t\sim N_a$ the autocorrelation turns slightly negative because of block-to-block anticorrelations in $\bar p_i$. Two characteristic times are then defined: the exponential correlation time $\tau_{O,a}$ from the exponential regime and a quasi-integrated time $\tau_{O,b}$ given by the maximum of the partial sum $\sum_{t'=0}^t\Gamma_O(t')$ before it crosses zero. Both are observed to scale as
$$
\tau_O\sim L^z
$$
[1004.1509].

With block size $N_a=100$ and width $\tilde v=0.1$, the measured dynamic exponents are, for $3$D Ising, $z_{m,a}=0.41(5)$, $z_{m,b}=0.38(3)$, $z_{e,a}=0.42(3)$, and $z_{e,b}=0.47(3)$; for $2$D $q=3$, they are $z_{m,a}=0.38(2)$, $z_{m,b}=0.29(2)$, $z_{e,a}=0.42(2)$, and $z_{e,b}=0.36(2)$. By contrast, the SW algorithm has $z\approx 0.46(3)$ in $3$D and $z\approx 0.49(1)$ in $2$D. The reported EIC exponents are also discussed relative to the Li–Sokal bound
$$
z\ge \alpha/\nu,
$$
with $\alpha/\nu\approx 0.172$ for $3$D Ising and $\alpha/\nu=2/5=0.4$ for $2$D $q=3$; within errors, the EIC values satisfy $z\gtrsim \alpha/\nu$, which is presented as confirmation that the dynamics remains an equilibrium cluster dynamics [1004.1509].

The auxiliary parameters are part of the algorithmic definition rather than mere implementation details. The block length $N_a$ must exceed the largest autocorrelation time on all $L$, roughly $N_a\gtrsim 2\,\max\tau_O$, so that block-to-block fluctuations in $\bar p_i$ remain subdominant. If $N_a$ is too small, one sees multiple sign-changes in $\Gamma_O(t)$ and a crossover to pure IC-like non-equilibrium behavior. The window width $v=\tilde v\,L^{-d/2}$ must match the equilibrium binomial variance. If $\tilde v$ is too large, the distribution of $p_{\{\sigma\}}$ broadens and the dynamics again crosses to IC-like behavior, with $z$ dropping well below $\alpha/\nu$; if $v$ decays faster than $L^{-d/2}$, the correlation time grows and $z$ increases [1004.1509].

## 5. Plaquette random-cluster generalization for Potts lattice gauge theory

A later development extends the invaded-cluster idea from the Potts spin model to the $q$-state Potts lattice-gauge theory on the $4$-torus $\mathbb T^4_N$. In this setting, the relevant objects are not FK bond-clusters but $2$-dimensional subcomplexes of plaquettes. The extension is carried out through a plaquette random-cluster model (PRCM), which plays for lattice gauge theory the role that the FK random-cluster model plays for the classical Potts model [2507.13503].

Let $X=\mathbb T^4_N$ be the $4$-dimensional cubical torus. The edge-spin configurations are
$$
C^1(X;\mathbb Z_q)=\{\text{edge-spin assignments } f:e\mapsto f(e)\in\mathbb Z_q\},
$$
and the plaquette coboundary is
$$
\delta^1f(x)=\sum_{e\subset\partial x}\pm\,f(e)\in\mathbb Z_q,
$$
for a plaquette $x$. The Potts lattice-gauge Gibbs measure at inverse temperature $\beta$ is
$$
\nu(f)=\frac{1}{Z(q,\beta)}\exp\Bigl(\beta\sum_{x\in X}1_{\{\delta^1f(x)=0\}}\Bigr).
$$
Introducing $p=1-e^{-\beta}$, one obtains an equivalent PRCM measure on subcomplexes
$$
P\subset X,\quad V_P=V_X,\;E_P=E_X,\;U_P\subset U_X
$$
given by
$$
\mu(P)\propto p^{\,|U_P|}(1-p)^{\,|U_X|-|U_P|}\,\bigl|H^1(P;\mathbb Z_q)\bigr|.
$$
Here $\bigl|H^1(P;\mathbb Z_q)\bigr|$ accounts for the gauge-invariance under $0$-cochains and replaces the usual $q^{\#\text{components}}$ factor in the FK random-cluster model [2507.13503].

This reformulation makes the generalization conceptually precise. In the spin-model IC algorithm, percolation is expressed by spanning $1$-cycles or wrapping clusters; in the lattice-gauge version, the corresponding critical event is homological percolation of $2$-cycles, that is, the appearance of spanning surfaces on the torus.

## 6. Homological stopping rule and empirical performance on the $4$-torus

In the generalized algorithm one maintains a joint state $(f_t,P_t)$, although only $f_t$ is passed from one Potts lattice-gauge update to the next. Given the current edge-spin configuration $f_t\in C^1(X;\mathbb Z_q)$, one generates a uniform random permutation of the plaquettes, initializes $P_0$ with all vertices and edges but no plaquettes, and then examines the plaquettes in that random order. A plaquette $x_k$ is added precisely when $\delta^1f_t(x_k)=0$; otherwise the current subcomplex is left unchanged. The process stops as soon as the invaded subcomplex contains a full set of nontrivial $2$-cycles in $H_2(\mathbb T^4_N;\mathbb F)$, after which one samples
$$
f_{t+1}\ \text{uniformly from}\ Z^1(P_{t+1};\mathbb Z_q)=\ker\bigl(\delta^1:C^1(P_{t+1};\mathbb Z_q)\to C^2(P_{t+1};\mathbb Z_q)\bigr)
$$
[2507.13503].

The stopping rule is expressed through persistent homology. If $\iota:P_k\hookrightarrow X$ denotes inclusion, the monitored group is
$$
PH_2(P_k)=\operatorname{im}\!\bigl(\iota_*:H_2(P_k;\mathbb F)\to H_2(X;\mathbb F)\bigr).
$$
On the $4$-torus, $H_2(\mathbb T^4_N;\mathbb F)$ has dimension $\binom{4}{2}=6$, and the algorithm stops when
$$
\operatorname{rank}PH_2(P_k)=6.
$$
Equivalently, the invaded subcomplex must contain six linearly independent spanning surfaces. This is the exact analogue, in degree $2$, of the winding-cluster criterion used in the classical IC method [2507.13503].

The critical point used in the simulations is the self-dual point
$$
p_c=\frac{\sqrt q}{1+\sqrt q},
$$
described as the conjectural critical point on $\mathbb T^4$. A theorem of Duncan–Schweinhart is also used to characterize the homological transition: if $P\sim \mu_{p,q}$ on $\mathbb T_N^{2d}$, and $R=\{\text{no } d\text{-cycles}\}$ and $S=\{\text{full basis of } d\text{-cycles}\}$, then as $N\to\infty$,
$$
\mu_{p,q}(R)\to 0\quad (p<p_c),\qquad \mu_{p,q}(S)\to 1\quad (p>p_c).
$$
This suggests that the homological stopping rule is sharply aligned with the critical regime in the thermodynamic limit.

The empirical data reported for $\mathbb T^4_{10}$ at $p=p_c$ compare Glauber dynamics, the plaquette-version SW algorithm, and the generalized invaded-cluster algorithm. For the energy
$$
\mathcal H(f_t)=-\sum_x1_{\{\delta^1f_t(x)=0\}},
$$
Glauber dynamics shows autocorrelation decay on a scale of thousands of iterations, whereas SW and invaded-cluster correlations fall to near zero in $\mathcal O(10)$ iterations. On a single high-end workstation, for $q=2$ with $N=10$ and $|U_X|=6000$ plaquettes, the per-iteration costs are approximately $0.57\,\mathrm s/\mathrm{it}$ for SW and approximately $1.7\,\mathrm s/\mathrm{it}$ for invaded-cluster using PHAT persistence over $\mathbb F_2$; for $q=3$, the corresponding figures are approximately $0.47\,\mathrm s/\mathrm{it}$ for SW and approximately $19\,\mathrm s/\mathrm{it}$ for invaded-cluster using Smith–normal-form persistence. Despite the higher per-iteration cost, invaded-cluster and SW mix in $\sim 10$ iterations versus $\sim 10^3$ for Glauber, giving a net speedup of $\sim 100\times$ in statistical efficiency for $\mathbb Z_2$ and similar for $\mathbb Z_3$. The simulations for $\mathbb Z_2$ and $\mathbb Z_3$ also indicate that the generalized algorithms allow efficient sampling on $4$-dimensional tori of linear scale at least $40$ [2507.13503].

The mechanism proposed for this acceleration is the same one that motivates cluster algorithms in the spin setting. Rather than flipping one edge at a time, both SW and invaded-cluster grow and then entirely re-sample an extended $2$-dimensional surface, or union of surfaces, whose geometry is dictated by the current frustration pattern $\{\delta^1f(x)=0\}$. This global move breaks up large connected regions of aligned plaquettes at criticality, dramatically reducing autocorrelation and circumventing the standard critical slowing-down of single-spin flips [2507.13503].

Source: https://www.emergentmind.com/topics/invaded-cluster-algorithm