Invaded-Cluster Algorithm in Statistical Physics

• Invaded-Cluster Algorithm is a Monte Carlo cluster method that grows clusters until a percolation criterion—such as a nontrivial winding or full basis of cycles—is met.
• The algorithm self-drives systems to criticality without needing a fixed bond probability, dramatically reducing critical slowing-down through nonlocal updates.
• Equilibriumlike variants and homological generalizations extend its application to canonical ensemble sampling and Potts lattice gauge theory on high-dimensional tori.

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 (Balog et al., 2010, Pizzimenti et al., 17 Jul 2025).

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) (Balog et al., 2010).

In SW dynamics one alternates a bond-placement step and a cluster-flip step. For each nearest-neighbor pair $(i,j)$ with $4$0, a bond is placed with probability $4$1; 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 $4$2, it erases all bonds, identifies the satisfied edges,

$4$3

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

$4$4

is then recorded, the bonds are erased, SW-clusters are rebuilt, and the clusters are flipped at random (Balog et al., 2010).

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 $4$5 converges to $4$6, the finite-size critical bond probability. Because each update is forced to stop at percolation, the method “self-drives” to the critical point (Balog et al., 2010).

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 $4$7 is far wider than the canonical binomial variance $4$8. As a consequence, the resulting ensemble is non-equilibrium and yields wrong finite-size scaling exponents (Balog et al., 2010). 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 $4$9 successive IC-style moves and constrains the fluctuations of $2$0 to a window of size $2$1 around the previous block’s mean $2$2, with

$2$3

At the end of block $2$4 one computes

$2$5

which becomes the center of the window for the next block (Balog et al., 2010).

The rationale is explicit. By choosing $2$6, one enforces

$2$7

matching the equilibrium binomial noise $2$8. The self-driving mechanism ensures $2$9, 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 (Balog et al., 2010).

The critical-scaling results reported for EIC reflect this restoration of equilibrium behavior. For the largest-cluster mass,

$q$0

the fits give $q$1 for the $q$2D Ising case with $q$3, and $q$4 for the $q$5D $q$6 Potts case with $q$7. For the thermal exponent obtained from

$q$8

the fits give $q$9 in $$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}.$$0D Ising and $$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}.$$1 in $$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}.$$2D $$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}.$$3 Potts. The corresponding three-parameter fits for the finite-size shift,

$$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}.$$4

yield $$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}.$$5 for $$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}.$$6D Ising and $$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}.$$7 for $$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}.$$8D $$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}.$$9. The reported exponents agree with known exact or best-numerical values to within $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$0-$$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$1 (Balog et al., 2010).

4. Dynamical behavior, autocorrelations, and parameter dependence

The dynamical analysis of EIC is formulated in terms of the normalized autocorrelation function

$$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$2

For $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$3, EIC exhibits a single-exponential decay,

$$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$4

while around $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$5 the autocorrelation turns slightly negative because of block-to-block anticorrelations in $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$6. Two characteristic times are then defined: the exponential correlation time $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$7 from the exponential regime and a quasi-integrated time $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$8 given by the maximum of the partial sum $$Z=\sum_{\gamma\in\Gamma}p^{b(\gamma)}(1-p)^{B-b(\gamma)}\,q^{c(\gamma)},$$9 before it crosses zero. Both are observed to scale as

$B$0

(Balog et al., 2010).

With block size $B$1 and width $B$2, the measured dynamic exponents are, for $B$3D Ising, $B$4, $B$5, $B$6, and $B$7; for $B$8D $B$9, they are $b(\gamma)$0, $b(\gamma)$1, $b(\gamma)$2, and $b(\gamma)$3. By contrast, the SW algorithm has $b(\gamma)$4 in $b(\gamma)$5D and $b(\gamma)$6 in $b(\gamma)$7D. The reported EIC exponents are also discussed relative to the Li–Sokal bound

$b(\gamma)$8

with $b(\gamma)$9 for $c(\gamma)$0D Ising and $c(\gamma)$1 for $c(\gamma)$2D $c(\gamma)$3; within errors, the EIC values satisfy $c(\gamma)$4, which is presented as confirmation that the dynamics remains an equilibrium cluster dynamics (Balog et al., 2010).

The auxiliary parameters are part of the algorithmic definition rather than mere implementation details. The block length $c(\gamma)$5 must exceed the largest autocorrelation time on all $c(\gamma)$6, roughly $c(\gamma)$7, so that block-to-block fluctuations in $c(\gamma)$8 remain subdominant. If $c(\gamma)$9 is too small, one sees multiple sign-changes in $(i,j)$0 and a crossover to pure IC-like non-equilibrium behavior. The window width $(i,j)$1 must match the equilibrium binomial variance. If $(i,j)$2 is too large, the distribution of $(i,j)$3 broadens and the dynamics again crosses to IC-like behavior, with $(i,j)$4 dropping well below $(i,j)$5; if $(i,j)$6 decays faster than $(i,j)$7, the correlation time grows and $(i,j)$8 increases (Balog et al., 2010).

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 $(i,j)$9-state Potts lattice-gauge theory on the $4$00-torus $4$01. In this setting, the relevant objects are not FK bond-clusters but $4$02-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 (Pizzimenti et al., 17 Jul 2025).

Let $4$03 be the $4$04-dimensional cubical torus. The edge-spin configurations are

$4$05

and the plaquette coboundary is

$4$06

for a plaquette $4$07. The Potts lattice-gauge Gibbs measure at inverse temperature $4$08 is

$4$09

Introducing $4$10, one obtains an equivalent PRCM measure on subcomplexes

$4$11

given by

$4$12

Here $4$13 accounts for the gauge-invariance under $4$14-cochains and replaces the usual $4$15 factor in the FK random-cluster model (Pizzimenti et al., 17 Jul 2025).

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

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

In the generalized algorithm one maintains a joint state $4$19, although only $4$20 is passed from one Potts lattice-gauge update to the next. Given the current edge-spin configuration $4$21, one generates a uniform random permutation of the plaquettes, initializes $4$22 with all vertices and edges but no plaquettes, and then examines the plaquettes in that random order. A plaquette $4$23 is added precisely when $4$24; otherwise the current subcomplex is left unchanged. The process stops as soon as the invaded subcomplex contains a full set of nontrivial $4$25-cycles in $4$26, after which one samples

$4$27

(Pizzimenti et al., 17 Jul 2025).

The stopping rule is expressed through persistent homology. If $4$28 denotes inclusion, the monitored group is

$4$29

On the $4$30-torus, $4$31 has dimension $4$32, and the algorithm stops when

$4$33

Equivalently, the invaded subcomplex must contain six linearly independent spanning surfaces. This is the exact analogue, in degree $4$34, of the winding-cluster criterion used in the classical IC method (Pizzimenti et al., 17 Jul 2025).

The critical point used in the simulations is the self-dual point

$4$35

described as the conjectural critical point on $4$36. A theorem of Duncan–Schweinhart is also used to characterize the homological transition: if $4$37 on $4$38, and $4$39 and $4$40, then as $4$41,

$4$42

This suggests that the homological stopping rule is sharply aligned with the critical regime in the thermodynamic limit.

The empirical data reported for $4$43 at $4$44 compare Glauber dynamics, the plaquette-version SW algorithm, and the generalized invaded-cluster algorithm. For the energy

$4$45

Glauber dynamics shows autocorrelation decay on a scale of thousands of iterations, whereas SW and invaded-cluster correlations fall to near zero in $4$46 iterations. On a single high-end workstation, for $4$47 with $4$48 and $4$49 plaquettes, the per-iteration costs are approximately $4$50 for SW and approximately $4$51 for invaded-cluster using PHAT persistence over $4$52; for $4$53, the corresponding figures are approximately $4$54 for SW and approximately $4$55 for invaded-cluster using Smith–normal-form persistence. Despite the higher per-iteration cost, invaded-cluster and SW mix in $4$56 iterations versus $4$57 for Glauber, giving a net speedup of $4$58 in statistical efficiency for $4$59 and similar for $4$60. The simulations for $4$61 and $4$62 also indicate that the generalized algorithms allow efficient sampling on $4$63-dimensional tori of linear scale at least $4$64 (Pizzimenti et al., 17 Jul 2025).

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 $4$65-dimensional surface, or union of surfaces, whose geometry is dictated by the current frustration pattern $4$66. 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 (Pizzimenti et al., 17 Jul 2025).