Worm Algorithm in Monte Carlo Methods

• Worm Algorithm is a Monte Carlo sampling technique that enlarges the state space by introducing a pair of defects to efficiently capture correlation functions and critical behaviors.
• It employs reversible and nonreversible update schemes across models like the Ising, gauge, and fermionic systems to optimize performance and reduce backtracking.
• The method provides rigorous mixing time bounds and adaptable estimators for nonlocal observables such as susceptibility, energy, and specific heat in diverse physical systems.

In the surveyed literature, the worm algorithm is presented as a versatile technique in the Markov chain Monte Carlo method for both classical and quantum systems, based on an enlarged state space with two defects, open boundaries, or open world-lines that are propagated locally and eventually removed or reclosed. Across Ising high-temperature expansions, dual flux or surface formulations, fermion-loop representations, and path-integral worldline methods, the defining operation is to move in an extended ensemble that relaxes a local closure constraint and then recover observables of the original closed sector from closed-worm visits or defect statistics [1703.03136] [2206.13881] [1409.4484].

1. Extended-ensemble construction

For the ferromagnetic Ising model on a bipartite lattice,
$$
H[\sigma] = -J\sum_{\langle ij\rangle}\sigma_i\sigma_j,\qquad \sigma_i=\pm1,\qquad K\equiv \beta J,
$$
the high-temperature expansion yields
$$
Z
= 2^N[\cosh K]^{N_b}
\sum_{{n_b\in{0,1}\ \text{forming only closed loops}}} (\tanh K)^\ell,
$$
with $\ell=\sum_b n_b$ and the even-degree constraint at each site. The worm construction enlarges this loop ensemble by temporarily allowing exactly two odd-degree sites, interpreted as defects or kinks. In the site-based formulation, the tail is at $x$, the head at $y$, and a head move $y\to y'$ flips one bond variable $n_b\to 1-n_b$ along $b=\langle y,y'\rangle$ [2206.13881].

On a general finite connected graph $G=(V,E)$, this structure is encoded by
$$
C_W := {A\subseteq E:\partial A=W},\qquad
W=C_0\cup C_2,
$$
where $\partial A$ is the set of odd-degree vertices of $(V,A)$. The worm state space is therefore the union of Eulerian subgraphs and subgraphs with exactly two odd vertices. In the zero-field ferromagnetic Ising model, the corresponding high-temperature weight is $\lambda(A)=x^{|A|}$ with $x=\tanh\beta$, and the worm chain is defined precisely on $C_0\cup C_2$ [1409.4484].

The same enlargement principle reappears outside spin models. In Abelian gauge theory, closed surfaces are supplemented by open-surface defects whose boundaries carry a Wilson-loop current [1011.1359]. In fermion-loop formulations of the $O(2N)$ Gross–Neveu model, one opens a fermionic string with two endpoints that sample a Majorana two-point function [1104.0569]. In canonical path-integral Monte Carlo for Bose systems, the diagonal sector of closed world-lines is enlarged by an off-diagonal sector in which one polymer is cut into a tail and a head [2203.00010]. The underlying mechanism is unchanged: a local constraint is violated in a controlled way, the defect pair explores configuration space, and closure restores the original sector.

2. Reversible worm formulations and observables

In the standard Ising worm update, a configuration with two kinks carries weight $\pi(C)\propto(\tanh K)^{\ell(C)}$. A head move $C\to C'$ changes $\ell\to \ell\pm1$, so the weight ratio is
$$
W(C\to C')\equiv \frac{\pi(C')}{\pi(C)}=\tanh K
\quad\text{or}\quad
\frac{1}{\tanh K}.
$$
With a symmetric proposal $T(C\to C')=1/(2d)$ to one of the $2d$ neighbors, the Metropolis acceptance
$$
A(C\to C')=\min{1,W(C\to C')}
$$
satisfies detailed balance. The basic life cycle is: insert two kinks at a random site, propagate the head by repeated flip-and-move steps, and remove the worm when the head returns to the tail [2206.13881].

In the lazy Metropolis formulation on a general graph, the chain again preserves detailed balance on $W=C_0\cup C_2$. One chooses a vertex $u$ uniformly in $V$ if the state lies in $C_0$, or uniformly in $\partial A$ if it lies in $C_2$, then chooses a neighboring vertex $v$, proposes $A'=A\Delta{uv}$, and accepts with a Metropolis–Hastings rule that includes the degree ratio $d(u)/d(v)$ when both current and proposed states lie in $C_2$ [1409.4484].

The reversible worm framework gives direct estimators for nonlocal observables. In the Ising chain on general graphs,
$$
\chi = \beta/\pi(C_0),\qquad
\langle \sigma_u\sigma_v\rangle = (n/2)\,\pi(C_{uv})/\pi(C_0)\quad (u\neq v),
$$
so susceptibility and two-point functions are encoded in the stationary weights of the closed and two-defect sectors [1409.4484]. In the two-dimensional Ising implementation used for graphene-impurity studies, the two-point correlator is obtained from the histogram of worm-end separations, while energy and specific heat are read from bond occupations [1303.0429]. The surveyed literature therefore treats the worm algorithm not as a direct spin-flip method, but as a sampler of graphical expansions whose defect statistics already carry the desired correlation information [1303.0429] [1703.03136].

3. Directed and lifted nonreversible variants

A central limitation of conventional worms is backscattering: local rejections or immediate reversals suppress diffusivity of the head and waste Monte Carlo time. Directed variants address this by attaching an incoming direction to the worm head and optimizing local scattering probabilities so that forward scattering is favored and backscattering is averted [1703.03136]. In the directed-worm algorithm described by Suwa, the kinks are placed on bonds rather than sites, the local bond state takes values $n_b\in{0,\tfrac12,1}$, and the kink carries a travel direction through the local scattering event [2206.13881].

The optimization is expressed as a geometric allocation problem. In the lifted setting one adds a discrete mode $\sigma\in{+1,-1}$ and defines stochastic flows
$$
v_{ij}^{\sigma\sigma'} \equiv \pi_i^\sigma\,p_{(i,\sigma)\to(j,\sigma')},
$$
subject to global balance
$$
\forall(i,\sigma):\qquad

\pi_i^\sigma

\sum_{j,\sigma'} v_{ij}^{\sigma\sigma'}

\sum_{j,\sigma'} v_{ji}^{\sigma'\sigma}.
$$
The allocation is then chosen to minimize backscattering flow, satisfy global balance, and maximize net stochastic flow from low-weight to high-weight scattering states. In the closed-form solution reported for the lifted directed-worm algorithm, one checks that $v_{\mathrm{back}}=0$ [2206.13881]. The same literature explicitly motivates lifting as a nonreversible perturbation that introduces net circulation in state space and cites mathematical proofs that any nonzero skew-detailed-balance perturbation shortens the relaxation time [2206.13881].

The numerical consequences are substantial. For the simple-cubic Ising model at the critical temperature, the directed worm is approximately $25$ times as efficient as the conventional worm update and more efficient than the Wolff cluster algorithm; the estimated dynamic critical exponent is $z\approx0.27$ [1703.03136]. For the four-dimensional hypercubic lattice Ising model at $1/T_c=0.1496947(5)$ and system sizes up to $L=56$, the lifted directed-worm algorithm yields
[
\text{asymp-var(Wolff)}/\text{asymp-var(present)}=5.4,\quad
\text{asymp-var(P--S worm)}/\text{asymp-var(present)}=80,\quad
\text{asymp-var(lifted B--S worm)}/\text{asymp-var(present)}=1.5
]
for the energy estimator, with corresponding factors $\simeq4.5$, $77$, and $1.7$ for $\chi$; the dynamic critical exponent for the energy is $z_{\mathrm{int}}^E\approx0$ [2206.13881]. Earlier lifted Ising worms also showed marked improvements: on toroidal grids in dimensions $d=2,\dots,5$, the asymptotic improvement factor $\tau_{\mathrm{int,rev}}/\tau_{\mathrm{int,irr}}$ grows from $2.4(4)$ in $d=2$ to $141(6)$ in $d=5$, while on the complete graph the lifted chain changes the scaling from $\tau_{\mathrm{int}}/n\approx\text{const}$ to $\tau_{\mathrm{int}}/n\sim n^{-1/2}$ [1711.05346].

4. Rigorous analysis and efficiency criteria

The reversible Prokof’ev–Svistunov worm algorithm for the zero-field ferromagnetic Ising model is not only practically effective but rigorously rapidly mixing on all finite connected graphs and at all temperatures. For $G=(V,E)$ with $|V|=n$, $|E|=m$, maximum degree $\Delta$, and $x=\tanh\beta\in(0,1)$, the lazy worm chain satisfies
$$
\tau_{\mathrm{mix}}(\delta)\le
\frac1{2x}\Bigl(\ln!\frac8x-\frac{\ln\delta}{m}\Bigr)
\Bigl(3+\frac1{mx}\Bigr)\Delta\,n^6\,m^2.
$$
In particular, $\tau_{\mathrm{mix}}=O(\Delta n^6 m^2\log n)$ is polynomial in graph size and $\log\delta^{-1}$. The same analysis yields fully polynomial randomized approximation schemes for the susceptibility and the two-point correlation function [1409.4484].

Performance comparisons in the worm literature are typically organized around autocorrelation and variance, rather than around acceptance rate alone. The standard definitions used in the directed-worm work include the exponential autocorrelation time $\tau_{\exp}$, the integrated autocorrelation time
$$
\tau_{\mathrm{int},\hat O}=\frac12+\sum_{t=1}^\infty A_{\hat O}(t),
$$
and the asymptotic variance
$$
v_{\mathrm{asymp},\hat O}=2\,\tau_{\mathrm{int},\hat O}\,\mathrm{Var}(\hat O),
$$
whose inverse is proportional to the sampling efficiency [1703.03136]. In the four-dimensional lifted directed-worm study, the normalized integrated autocorrelation time is measured in units of $N=L^4$ spin-updates, the estimator variance is written as $v=\mathrm{Var}[\hat O]/\langle\hat O\rangle^2$, and the asymptotic variance is reported as $v_{\mathrm{asymp}}=2\tau_{\mathrm{int}}v$ [2206.13881].

This metric choice matters because different observables can display different dynamic behavior. In the three-dimensional directed-worm study, the worm and Wolff algorithms produce different exponents for the integrated autocorrelation time of the magnetic susceptibility estimator but the same exponent for the asymptotic variance [1703.03136]. The literature therefore treats efficiency as an observable-dependent quantity that combines decorrelation and estimator variance.

5. Generalizations across models and representations

In Abelian gauge theory, the worm idea extends from closed loops to closed surfaces. For the lattice $U(1)$ gauge theory with Wilson plaquette action, one samples the all-order strong-coupling expansion in integer-valued plaquette occupation numbers and enlarges the ensemble to include open-surface defects whose boundaries are Wilson loops. In three dimensions at fixed physical volume, the worm-inspired algorithm exhibits $\tau_{\mathrm{int}}\approx O(1)$ essentially independent of $L$, whereas the Metropolis $\tau_{\mathrm{int}}$ grows roughly $\propto L^z$ with $z\approx2$; the dimensionless cost indicator is reduced by a factor $3$–$5$ for the average plaquette and by an order of magnitude for the effective mass at $L=40$ [1011.1359]. The surface worm algorithm for Abelian Gauge–Higgs systems generalizes the original loop worm to dual configurations of surfaces and loops, outperforms a local update for a wide range of parameters, and in the $Z(3)$ case at finite chemical potential avoids the complex-action problem because the dual weights remain positive [1211.3436] [1304.3278].

For scalar and sigma-like systems, worm methods remain effective outside the strict Ising limit. In one-component lattice $\phi^4$ theory, the link-variable worm algorithm works from $\lambda=0$ to $\lambda=\infty$; at $\lambda=0$ and fixed $mL=4$ it shows moderate critical slowing down, with, for example, $\tau_{\mathrm{int}}(E)\propto L^{1.62(2)}$ and $\tau_{\mathrm{int}}(m)\propto L^{1.86(7)}$ in $D=2$, while at $\lambda=\tfrac12$ the growth is very mild and close to $A+B\ln L$, and at $\lambda=\infty$ earlier studies show virtually no critical slowing down [1101.3452]. In the complex $\phi^4$ model, the worm framework can sample not only the charged correlator $\langle\phi^*(x)\phi(y)\rangle$, but also $\langle|\phi(x)||\phi(y)|\rangle$, $\langle\arg(\phi(x))\arg(\phi(y))\rangle$, and condensates such as $\langle|\phi|\rangle$ at every step of the Monte Carlo evolution rather than only on closed-worm configurations [1602.09017]. For the two-dimensional $\mathrm{CP}(N-1)$ model, dual flux-variable worms and internal-space sub-worm moves permit simulations at finite density without sign problems; the internal-space moves restore ergodicity for $N>2$, and the quartic and auxiliary-$U(1)$ lattice actions show marked differences in their approach to the continuum limit [1610.01435] [1703.08571].

Loop, fermionic, and quantum worldline formulations admit equally direct worm constructions. A worm algorithm for the $O(n)$ loop model on any finite connected bipartite cubic graph works for any real $n>0$ and any edge weight, including the fully-packed limit, and is rigorously ergodic with the correct stationary distribution; at $n=2$, it yields simulations of zero-temperature Potts antiferromagnets for which the Wang–Swendsen–Kotecký algorithm is non-ergodic [1011.1980]. In the $O(2N)$ Gross–Neveu model with Wilson fermions, the open fermionic string samples two-point functions directly, extends to bound states involving an arbitrary number of Majorana fermions, and eliminates critical slowing down in practice [1104.0569]. In canonical path-integral Monte Carlo for Bose systems with periodic boundary conditions, the worm algorithm is formulated so that only the initial point of each path lies inside the simulation box while the remaining points span the whole space, allowing detailed balance without ambiguity in the periodic image selection; benchmarks against exact finite-$N$ results for the ideal Bose gas confirm the construction [2203.00010]. In the quantum transverse-field Ising model, a continuous-time worm algorithm in the field-direction basis yields the two-dimensional critical point $h_c=3.044\,330(6)$ and loop fractal dimensions $d_{\ell\downarrow}(1\mathrm{D})=1.37(1)\approx11/8$ and $d_{\ell\downarrow}(2\mathrm{D})=1.75(3)$ [2005.10066].

6. Scope, misconceptions, and recent directions

The surveyed literature indicates that “worm algorithm” denotes a family of constructions rather than one unique transition rule. Reversible Metropolis worms, directed worms, lifted worms, surface worms, internal-space sub-worms, open-string fermion worms, and path-integral worms all share the defect-enlargement strategy, but they differ in local state spaces, balance conditions, and observables. Another recurrent misconception is that worms are intrinsically tied to detailed balance: the standard Prokof’ev–Svistunov form is reversible, whereas lifting and directed scattering deliberately break detailed balance to introduce net stochastic flow and suppress backtracking [2206.13881] [1711.05346].

A further misconception is that worms are useful only for two-point functions. The corpus surveyed here uses worm statistics to estimate susceptibility, energy, specific heat, Wilson loops, plaquette correlators, effective masses, condensates, and more general composite correlators [1011.1359] [1602.09017]. This breadth follows directly from the enlarged ensemble: once the defect sector is built to represent an insertion of the relevant operator, the same Markov chain can update configurations and sample the correlator.

A recent cross-disciplinary development applies the same $C_0\cup C_2$ logic to decoding. For matchable qLDPC codes, the worm algorithm approximately computes probabilities of logical error classes given a syndrome, thereby implementing approximate optimal decoding. The underlying Markov chain has a rigorous relaxation-time bound in terms of the defect susceptibilities $\chi_2=Z_2/Z_0$ and $\chi_4=Z_4/Z_0$, and the method applies to the surface code, the honeycomb Floquet code, and hyperbolic surface codes with and without measurement errors. Numerically, the worm decoder gives $p_c^{\rm worm}=0.0310\pm0.0005$ for the rotated surface code with independent bit-flip and measurement error, $p_c^{\rm worm}=0.019\pm0.001$ for a family of hyperbolic ${5,5}$ codes, and a correlated-worm threshold of approximately $0.186$ for depolarizing noise on the rotated surface code, compared to approximately $0.154$ for plain minimum-weight perfect matching and approximately $0.160$ for correlated matching [2603.05428]. This suggests that the worm paradigm has evolved from a specialized update for lattice loop gases into a general framework for sampling constrained subgraph distributions.