% % Created by Jaron Kent-Dobias on Thu Apr 20 12:50:56 EDT 2017. % Copyright (c) 2017 Jaron Kent-Dobias. All rights reserved. % \documentclass[aps,prl,reprint]{revtex4-1} \usepackage[utf8]{inputenc} \usepackage{amsmath,amssymb,latexsym,mathtools,xifthen,upgreek} \usepackage{algcompatible} \usepackage{algorithm} % uncomment to label only equations that are referenced in the text %\mathtoolsset{showonlyrefs=true} % I want labels but don't want to type out ``equation'' \def\[{\begin{equation}} \def\]{\end{equation}} % math not built-in \def\arcsinh{\mathop{\mathrm{arcsinh}}\nolimits} \def\arccosh{\mathop{\mathrm{arccosh}}\nolimits} \def\ei{\mathop{\mathrm{Ei}}\nolimits} % exponential integral Ei \def\re{\mathop{\mathrm{Re}}\nolimits} \def\im{\mathop{\mathrm{Im}}\nolimits} \def\tr{\mathop{\mathrm{Tr}}\nolimits} \def\sgn{\mathop{\mathrm{sgn}}\nolimits} \def\dd{d} % differential \def\O{O} % big O \def\o{o} % little O \def\Z{\mathbb Z} \def\R{\mathbb R} \def\G{S} \def\X{\mathbb X} \def\inf{\mathrm{inf}} % subscript for ``critical'' values, e.g., T_\c \def\c{\mathrm c} % scaling functions \def\fM{\mathcal M} % magnetization \def\fX{\mathcal Y} % susceptibility \def\fF{\mathcal F} % free energy \def\fiF{\mathcal H} % imaginary free energy \def\fS{\mathcal S} % surface tension \def\fG{\mathcal G} % exponential factor \def\H{\mathcal H} % lattice types \def\sq{\mathrm{sq}} \def\tri{\mathrm{tri}} \def\hex{\mathrm{hex}} % dimensions \def\dim{d} \def\twodee{\textsc{2d} } \def\threedee{\textsc{3d} } \def\fourdee{\textsc{4d} } % fancy partial derivative \newcommand\pd[3][]{ \ifthenelse{\isempty{#1}} {\def\tmp{}} {\def\tmp{^#1}} \frac{\partial\tmp#2}{\partial#3\tmp} } \newcommand\nn[2]{\langle#1#2\rangle} \newcommand\avg[1]{\langle#1\rangle} \newcommand\eavg[1]{\langle#1\rangle_{\mathrm e}} \newcommand\mavg[1]{\langle#1\rangle_{\mathrm m}} \def\e{\mathrm e} \def\m{\mathrm m} \newcommand\new[1]{\tilde{#1}} \renewcommand\vec[1]{\mathbf{#1}} \newcommand\unit[1]{\hat{\mathbf{#1}}} \def\mcmc{\textsc{McMC}} \def\vsigma{\pmb{\upsigma}} % used to reformat display math to fit in two-column ``reprint'' mode \makeatletter \newif\ifreprint \@ifclasswith{revtex4-1}{reprint}{\reprinttrue}{\reprintfalse} \makeatother \begin{document} \title{An efficient cluster algorithm for spin systems in a symmetry-breaking field} \author{Jaron Kent-Dobias} \author{James P.~Sethna} \affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA} \date\today \begin{abstract} We introduce a generalization of the `ghost spin' representation of spin systems that restores full symmetry group invariance in an arbitrary external field via the introduction of a `ghost transformation.' This offers a natural way to extend celebrated spin-cluster Monte Carlo algorithms to systems in arbitrary fields by running the ordinary cluster-flipping process on the new representation. For several canonical systems, we show that this extension with field preserves the scaling of dynamics so celebrated without field. \end{abstract} \maketitle Spin systems are important in the study of statistical physics and phase transitions. Rarely exactly solvable, they are typically studied by approximation methods and numeric means. Monte Carlo methods are a common way of doing this, approximating thermodynamic quantities by sampling the distribution of systems states. For a particular system, a Monte Carlo algorithm is better the faster it arrives at a statistically independent sample. This is typically a problem at critical points, where critical slowing down \cite{wolff_critical_1990} results in power-law divergences of any dynamics. Celebrated cluster algorithms largely addressed this for many spin systems in the absence of external fields by using nonlocal updates \cite{janke_nonlocal_1998} whose clusters undergo a percolation transition at the critical point of the system \cite{coniglio_clusters_1980} and that in relatively small dynamic exponents \cite{wolff_comparison_1989,du_dynamic_2006,liu_dynamic_2014,wang_cluster_1990}, including the Ising, $\mathrm O(n)$ \cite{wolff_collective_1989}, and Potts \cite{swendsen_nonuniversal_1987,baillie_comparison_1991} models. These algorithms rely on the natural symmetry of the systems in question under global rotations, so the general application of external fields is not trivial. Some success has been made in extending these algorithms to systems in certain external fields based on applying the ghost site representation \cite{coniglio_exact_1989} of certain spin systems that returns global rotation invariance to spin Hamiltonians at the cost of an extra degree of freedom, but these results only allow the application of a narrow category of fields \cite{alexandrowicz_swendsen-wang_1989,destri_swendsen-wang_1992,lauwers_critical_1989,wang_clusters_1989}. We show that the scaling of correlation time near the critical point of several models suggests that this approach is a natural one, e.g., that it extends the celebrated scaling of dynamics in these algorithms at zero field to various non-symmetric perturbations. We also show, by a redefinition of the spin--spin coupling in a generic class of such systems, systems with arbitrary external fields applied can be treated using cluster methods. Let $G=(V,E)$ be a graph, where the set of vertices $V=\{1,\ldots,N\}$ enumerates the sites of a lattice and the set of edges $E$ contains pairs of neighboring sites. Let $R$ be a group acting on a set $X$, with the action of group elements $r\in R$ on elements $s\in X$ denoted $r\cdot s$. $X$ is the set of states accessible by a spin, and $R$ is the \emph{symmetry group} of $X$. The set $X$ must admit a measure $\mu$ that is invariant under the action of $R$, e.g., for any $A\subseteq X$ and $r\in R$, $\mu(r\cdot A)=\mu(A)$. This trait is shared by the counting measure on any discrete set, or by any group acting by isometries on a Riemannian manifold, such as $\mathrm O(n)$ on $S^{n-1}$ in the $\mathrm O(n)$ model \cite{caracciolo_wolff-type_1993}. Finally, the subset of elements in $R$ of order two must act transitively on $X$. This property, while apparently obscure, is shared by any symmetric space \cite{loos_symmetric_1969} or by any transitive, finitely generated isometry group. In fact, all the examples listed here have spins spaces with natural metrics whose symmetry group is the set of isometries of the spin spaces. We put one spin at each site of the lattice described by $G$, so that the state of the entire spin system is described by elements $\vec s\in X\times\cdots\times X=X^N$. The Hamiltonian of this system is a function $\H:X^N\to\R$ defined by \[ \H(\vec s)=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!Z(s_i,s_j)-\sum_{i\in V}B(s_i), \] where $Z:X\times X\to\R$ couples adjacent spins and $B:X\to\R$ is an external field. $Z$ must be symmetric in its arguments and invariant under the action of any element of $R$ applied to the entire lattice, that is, for any $r\in R$ and $s,t\in X$, $Z(r\cdot s,r\cdot t)=Z(s,t)$. One may also allow $Z$ to also be a function of the edge---for modelling random-bond, long-range, or anisotropic interactions---or allow $B$ to be a function of site---for applying arbitrary boundary conditions or modelling random fields. All the formal results of this paper hold equally well for these cases, but we will drop the additional index notation for clarity. \begin{table*}[htpb] \begin{tabular}{l||ccccc} & Spins ($X$) & Symmetry ($R$) & Action ($g\cdot s$) & Coupling ($Z(s,t)$) & Common Field ($B(s)$) \\ \hline\hline Ising & $\{-1,1\}$ & $\Z/2\Z$ & $0\cdot s\mapsto s$, $1\cdot s\mapsto -s$ & $st$ & $Hs$ \\ $\mathrm O(n)$ & $S^{n-1}$ & $\mathrm O(n)$ & $M\cdot s\mapsto Ms$ & $s^{\mathrm T}t$ & $H^{\mathrm T}s$\\ Potts & $\mathbb Z/q\mathbb Z$ & $D_n$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $\delta(s,t)$ & $\sum_mH_m\delta(m,s)$\\ Clock & $\mathbb Z/q\mathbb Z$ & $D_n$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $\cos(2\pi\frac{s-t}q)$ & $\sum_mH_m\cos(2\pi\frac{s-m}q)$\\ Discrete Gaussian & $\mathbb Z$ & $D_\inf$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $(s-t)^2$ & $Hs^2$\\ \end{tabular} \caption{Several examples of spin systems and the symmetry groups that act on them. Common choices for the spin--spin coupling in these systems and their external fields are also given. Other fields are possible, of course: for instance, some are interested in modulated fields $H\cos(2\pi k\theta(s))$ for integer $k$ and $\theta(s)$ giving the angle of $s$ to some axis applied to $\mathrm O(n)$ models \cite{jose_renormalization_1977}.} \label{table:models} \end{table*} The goal of statistical mechanics as applied to these systems is to compute expectation values of observables $A:X^N\to\R$. Assuming the ergodic hypothesis holds (for systems with broken-symmetry states, it does not), the expected value $\avg A$ of an observable $A$ is its average over every state $\vec s$ in the configuration space $X^N$ weighted by the probability $p(\vec s)$ of that state appearing, or \[ \avg A =\frac{\int_{X^N}A(\vec s)e^{-\beta\H(\vec s)}\,\dd\mu(\vec s)} {\int_{X^N}e^{-\beta\H(\vec s)}\,\dd\mu(\vec s)} \] where for $Y_1\times\cdots\times Y_N=Y\subseteq X^N$ the measure $\mu(Y)=\mu(Y_1)\cdots\mu(Y_N)$ is the simple extension of the measure on $X$ to a measure on $X^N$. These values are estimated by Monte Carlo techniques by constructing a finite sequence of states $\{\vec s_1,\ldots,\vec s_M\}$ such that \[ \avg A\simeq\frac1M\sum_{i=1}^MA(\vec s_i) \] Sufficient conditions for this average to converge to $\avg A$ as $M\to\infty$ are that the process that selects $\vec s_{i+1}$ given the previous states be Markovian (only depends on $\vec s_i$), ergodic (any state can be accessed), and obey detailed balance (the ratio of probabilities that $\vec s'$ follows $\vec s$ and vice versa is equal to the ratio of weights for $\vec s$ and $\vec s'$ in the ensemble). While any several related cluster algorithms can be described for this system, we will focus on the Wolff algorithm in particular \cite{wolff_collective_1989}. We will first describe a generalized version of the celebrated Wolff algorithm in the standard case where $B(s)=0$. After reflecting on the technical requirements of that algorithm, we will introduce a transformation to our system and Hamiltonian that allows the same algorithm to be applied with nonzero, in fact \emph{arbitrary}, external fields. The Wolff algorithm proceeds in the following way. \begin{enumerate} \item Pick a random site and a random rotation $r\in R$ of order two, and add the site to a stack. \item While the stack isn't empty, \begin{enumerate} \item pop site $m$ from the stack. \item If site $m$ isn't marked, \begin{enumerate} \item mark the site. \item For every $j$ such that $\{m,j\}\in E$, add site $j$ to the stack with probability \[ p_r(s_m,s_j)=\min\{0,1-e^{\beta(Z(r\cdot s_m,s_j)-Z(s_m,s_j))}\}. \] \item Take $s_m\mapsto r\cdot s_m$. \end{enumerate} \end{enumerate} \end{enumerate} When the stack is exhausted, a cluster of connected spins will have been rotated by the action of $r$. In order for this algorithm to be useful, it must satisfy ergodicity and detailed balance. The probability $P(\vec s\to\vec s')$ that the configuration $\vec s$ is brought to $\vec s'$ by the flipping of a cluster formed by accepting rotations of spins via bonds $C\subseteq E$ and rejecting rotations via bonds $\partial C\subset E$ is related to the probability of the reverse process $P(\vec s'\to\vec s)$ by \begin{widetext} \[ \begin{aligned} \frac{P(\vec s\to\vec s')}{P(\vec s'\to\vec s)} &=\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r^{-1}}(s_i',s_j')}\prod_{\{i,j\}\in\partial C}\frac{1-p_r(s_i,s_j)}{1-p_{r^{-1}}(s'_i,s'_j)} =\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r}(r\cdot s_i,r\cdot s_j)}\prod_{\{i,j\}\in\partial C}\frac{1-p_r(s_i,s_j)}{1-p_{r^{-1}}(r\cdot s_i,s_j)} \\ &=\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r}(s_i,s_j)}\prod_{\{i,j\}\in\partial C}e^{\beta(Z(r\cdot s_i,s_j)-Z(s_i,s_j))} =\frac{e^{-\beta\H(\vec s)}}{e^{-\beta\H(\vec s')}} \end{aligned} \] \end{widetext} whence detailed balance is satisfied. Ergodicity is satisfied since we have ensured that the subset of elements in $R$ that are order two acts transitively on $K$, e.g., for any $s,t\in X$ there exists $r\in R$ such that $r\cdot s=t$. Since there is a nonzero probability that only one spin is rotated and that spin can be rotated into any state, ergodicity follows. The function of the algorithm described above depends on the fact that the coupling $Z$ depends only on the relative orientation of the spins---global reorientations by acting by some rotation do not affect the Hamiltonian. The external field $B$ breaks this symmetry. However, this can be resolved. Define a new graph $\tilde G=(\tilde V,\tilde E)$, where $\tilde V=\{0,1,\ldots,N\}$ adds a new `ghost' site $0$ which is connected by \[ \tilde E=E\cup\big\{\{0,i\}\mid i\in V\big\} \] to all other sites. Instead of assigning this ghost site a spin whose value comes from the set $X$, we will assign it values in the symmetry group $s_0\in R$, so that the new configuration space of the model is $R\times X^N$. We introduce a Hamiltonian $\tilde\H:R\times X^N\to\R$ defined by \[ \begin{aligned} \tilde\H(s_0,\vec s) &=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!Z(s_i,s_j) -\sum_{i\in V}B(s_0^{-1}\cdot s_i)\\ &=-\!\!\!\!\sum_{\{i,j\}\in\tilde E}\!\!\!\!\tilde Z(s_i,s_j) \end{aligned} \] where the new coupling $\tilde Z:(R\cup X)\times(R\cup X)\to\R$ is defined for $s,t\in R\cup X$ by \[ \tilde Z(s,t) = \begin{cases} Z(s,t) & \text{if $s,t\in X$} \\ B(s^{-1}\cdot t) & \text{if $s\in R$} \\ B(t^{-1}\cdot s) & \text{if $t\in R$} \end{cases} \label{eq:new.z} \] Note that this modified coupling is invariant under the action of group elements: for any $r,s_0\in R$ and $s\in X$, \[ \begin{aligned} \tilde Z(rs_0,r\cdot s) &=B((rs_0)^{-1}\cdot (r\cdot s))\\ &=B((s_0^{-1}r^{-1})\cdot(r\cdot s))\\ &=B((s_0^{-1}r^{-1}r)\cdot s)\\ &=B(s_0^{-1}\cdot s) =\tilde Z(s_0,s) \end{aligned} \] The invariance $\tilde Z$ to rotations given other arguments follows from the invariance properties of $Z$. We have produced a system that incorporates the field function $B$ whose Hamiltonian is invariant to global rotations, but how does it relate to our previous system, whose properties we actually want to measure? If $A:X^N\to\R$ is an observable of the original system, one can construct an observable $\tilde A:R\times X^N\to\R$ of the new system defined by \[ \tilde A(s_0,\vec s)=A(s_0^{-1}\cdot\vec s) \] whose expectation value in the new system equals that of the original observable in the old system. First, note that $\tilde\H(1,\vec s)=\H(\vec s)$. Since the Hamiltonian is invarient under global rotations, it follows that for any $g\in R$, $\tilde\H(g,g\cdot\vec s)=\tilde\H(g^{-1}g,g^{-1}g\cdot\vec s)=\tilde\H(1,\vec s)=\H(\vec s)$. Using the invariance properties of the measure on $X$ and introducing a measure $\rho$ on $R$, it follows that \[ \begin{aligned} \avg{\tilde A} &=\frac{ \int_R\int_{X^N}\tilde A(s_0,\vec s)e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(s_0^{-1}\cdot\vec s)e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(\vec s')e^{-\beta\tilde\H(s_0,s_0\cdot\vec s')}\dd\mu(s_0\cdot\vec s')\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,s_0\cdot\vec s')}\dd\mu(s_0\cdot\vec s')\,\dd\rho(s_0) }\\ &=\frac{ \int_R\dd\rho(s_0)}{ \int_R\dd\rho(s_0)}\frac{\int_{X^N}A(\vec s')e^{-\beta\H(\vec s')}\dd\mu(\vec s') }{\int_{X^N}e^{-\beta\H(\vec s')}\dd\mu(\vec s') } =\avg A \end{aligned} \] To summarize, spin systems in a field may be treated in the following way. \begin{enumerate} \item Add a site to your lattice adjacent to every other site. \item Initialize a ``spin'' at that site that is a representation of a member of the symmetry group of your ordinary spins. \item Carry out the ordinary Wolff cluster-flip procedure on this new lattice, substituting $\tilde Z$ as defined in \eqref{eq:new.z} for $Z$. \end{enumerate} Ensemble averages of observables $A$ can then be estimated by sampling the value of $\tilde A$ on the new system. In contrast with the simpler ghost spin representation, this form of the Hamiltonian mya be considered the ``ghost transformation'' representation. \section{Examples} \subsection{The Ising Model} In the Ising model, spins are drawn from the set $\{1,-1\}$. The symmetry group of this model is $C_2$, the cyclic group on two elements, which can be conveniently represented by the multiplicative group with elements $\{1,-1\}$, exactly the same as the spins themselves. The only nontrivial element is of order two. Because the symmetry group and the spins are described by the same elements, performing the algorithm on the Ising model in a field is very accurately described by simply adding an extra spin coupled to all others and running the ordinary algorithm. The ghost spin version of the algorithm has been applied by several researchers previously \cite{wang_clusters_1989,ray_metastability_1990,destri_swendsen-wang_1992,lauwers_critical_1989} \subsection{The $\mathrm O(n)$ Model} In the $\mathrm O(n)$ model, spins are described by vectors on the $(n-1)$-sphere, so that $X=S^{n-1}$. The symmetry group of this model is $O(n)$, $n\times n$ orthogonal matrices. The symmetry group acts on the spins by matrix multiplication. The elements of $O(n)$ that are order two are reflections about some hyperplane through the origin and $\pi$ rotations about any axis through the origin. Since the former generate the entire group, the set of reflections alone suffices to provide ergodicity. Computation of the coupling of ordinary spins with the external field and expectation values requires a matrix inversion, but since the matrices in question are orthogonal this is quickly accomplished by a transpose. The ghost-spin version of the algorithm has been used to apply a simple vector field by previous researchers \cite{dimitrovic_finite-size_1991}. \subsection{The Potts \& Clock Models} In both the $q$-state Potts and clock models, spins are described by $\Z/q\Z$, the set of integers modulo $q$. The symmetry group of this model is the dihedral group $D_q=\{r_0,\ldots,r_{q-1},s_0,\ldots,s_{q-1}\}$, the group of symmetries of a regular $q$-gon. The element $r_n$ represents a rotation of the polygon by $2\pi n/q$, and the element $s_n$ represents a reflection composed with a rotation $r_n$. The group acts on the spins by permutation: $r_n\cdot m={n+m}\pmod q$ and $s_n\cdot m={-(n+m)}\pmod q$. Intuitively, this can be thought of as the natural action of the group on the vertices of a regular polygon that have been numbered $0$ through $q-1$. The elements of $D_q$ that are of order 2 are all reflections and $r_{q/2}$ if $q$ is even, though the former can generate the latter. While the reflections do not necessarily generate the entire group, for any $n,m\in\Z/q\Z$ there exists a reflection that takes $n\to m$, ensuring ergodicity. The elements of the dihedral group can be stored simply as an integer and a boolean that represents whether the element is a pure rotation or a reflection. The principle difference between the Potts and clock models is that, in the latter case, the form of the coupling $Z$ allows a geometric interpretation as being two-dimensional vectors fixed with even spacing along the unit circle. \subsection{Discrete (or Continuous) Gaussian Model} Though not often thought of as a spin model, simple roughening of surfaces can be described in this framework. The set of states is the integers $\Z$ and its symmetry group is the infinite dihedral group $D_\infty=\{r_i,s_i\mid i\in\Z\}$, where the action of the symmetry on the spins $j\in\Z$ is given by $r_i\cdot j=i+j$ and $s_i\cdot j=-i-j$. These are shifts by $i$ and reflection about the integer $i$, respectively. The elements of order two are the reflections $s_i$, which suffice to provide ergodicity as any integer can be taken to any other in one step of this kind. The coupling is usually taken to be $Z(i,j)=(i-j)^2$, though it may also be any function of the absolute difference $|i-j|$. Because random choices of integer will almost always result in energy changes so big that the whole system is always flipped, it is better to select random reflections about integers close to the average state of the system. Continuous roughening models---where the spin states are described by real numbers and the symmetry group is $\mathrm E(1)$, the Euclidean group for one-dimensional space---are equally well described. A variant of the algorithm has been applied without a field before \cite{evertz_stochastic_1991}. %\begin{figure} % \centering % \input{fig_correlation} % \caption{The autocorrelation time $\tau$ of the internal energy $\H$ for a $n=32\times32$ square-lattice Ising % model simulated by the three nonzero field methods detailed in this paper. % The top plot is for simulations at high temperature, $T=5.80225$, the middle %plot is for simulations at the critical temperature, $T=2.26919$, and the %bottom plot is for simulations at } % \label{fig:correlation} %\end{figure} \section{Dynamic scaling} No algorithm worthwhile if it doesn't run efficiently. Our algorithm, being an extension of the Wolff algorithm, should be considered successful if it likewise extends its efficiency in the systems that algorithm succeeds. The Wolff algorithm succeeds at Cluster algorithms were celebrated for their small dynamic exponents $z$, which with the correlation time $\tau$ scales like $L^z$, where $L=N^{-D}$. In the vicinity of the critical point, the renormalization group predicts scaling behavior for the correlation time of the form \[ \tau=t^{-z\nu}\mathcal T(ht^{-\beta\delta},Lt^\nu) =h^{-z\nu/\beta\delta}\mathcal T'(ht^{-\beta\delta},Lt^\nu). \] If a given dynamics for a system at zero field results in scaling like $t^{-z\nu}$, one should expect its natural extension in the presence of a field to scale like $h^{-z\nu/\beta\delta}$. We measured the autocorrelation time for the 2D square-lattice model at a variety of system sizes, temperatures, and fields $B(s)=hs/\beta$ using methods here \cite{geyer_practical_1992}. The resulting scaling behavior, plotted in Fig.~\ref{fig:correlation_time-collapse}, is indeed consistent with the zero-field scaling behavior. \begin{figure} \centering \input{fig_correlation_collapse-hL} \input{fig_correlation-temp} \caption{Collapses of the correlation time $\tau$ of the 2D square lattice Ising model (top) along the critical isotherm at various systems sizes $N=L\times L$ for $L=8$, $16$, $32$, $64$, $128$, and $256$ as a function of the renormalization invariant $hL^{\beta\delta/\nu}$ and (bottom) in the low-temperature phase at $L=128$ for various temperatures as a function of the invariant $ht^{-\beta\delta}$. } \label{fig:correlation_time-collapse} \end{figure} Since the formation and flipping of clusters is the hallmark of the Wolff dynamics, another way to ensure that the dynamics with field scale like those without is to analyze the distribution of cluster sizes. The success of the algorithm at zero field is related to the way that clusters formed undergo a percolation transition at models' critical point. According to the scaling theory of percolation \cite{stauffer_scaling_1979}, the distribution of cluster sizes in a full decomposition of the system scales consistently near the critical point if it has the form \[ P_{\text{SW}}(s)=s^{-\tau}f(ts^\sigma,th^{-1/\beta\delta},tL^{1/\nu}). \] The distribution of cluster sizes in the Wolff algorithm can be computed from this using the fact that the algorithm selects clusters with probability proportional to their size, or \[ \begin{aligned} \avg{s_{\text{\sc 1c}}}&=\sum_ssP_{\text{\sc 1c}}(s)=\sum_ss\frac sNP_{\text{SW}}(s)\\ &=t^{-\gamma}g(th^{-1/\beta\delta},tL^{1/\nu})\\ &=L^{\gamma/\nu}\mathcal G(ht^{-\beta\delta},hL^{\beta\delta/\nu}) \end{aligned} \] For the Ising model, an additional scaling relation can be written. Since in that case the average cluster size is the average squared magnetization, it can be related to the scaling functions of the magnetization and susceptibility per site by (with $ht^{-\beta\delta}$ dependence dropped) \[ \begin{aligned} \avg{s_{\text{\sc 1c}}} &=L^{D}\avg{M^2}=\beta\avg\chi+L^{D}\avg{M}^2\\ &=L^{\gamma/\nu}\big[(hL^{\beta\delta/\nu})^{-\gamma/\beta\delta}\beta \mathcal Y(hL^{\beta\delta/\nu})\\ &\hspace{7em}+(hL^{\beta\delta/\nu})^{2/\delta}\mathcal M(hL^{\beta\delta/\nu})\big]. \end{aligned} \] We therefore expect that, for the Ising model, $\avg{s_{\text{\sc 1c}}}$ should go as $(hL^{\beta\delta})^{2/\delta}$ for large argument. We further conjecture that this scaling behavior should hold for other models whose critical points correspond with the percolation transition of Wolff clusters. This behavior is supported by our numeric work along the critical isotherm for various Ising, Potts, and $\mathrm O(n)$ models, shown in Fig.~\ref{fig:cluster_scaling}. Fields for the Potts and $\mathrm O(n)$ models take the form $B(s)=(h/\beta)\sum_m\cos(2\pi(s-m)/q)$ and $B(s)=(h/\beta)[1,0,\ldots,0]s$ respectively. As can be seen, the average cluster size collapses for each model according to the scaling hypothesis, and the large-field behavior likewise scales as we expect from the na\"ive Ising conjecture. \begin{figure*} \input{fig_clusters_ising2d} \caption{Collapses of rescaled average Wolff cluster size $\avg s_{\text{\sc 1c}}L^{-\gamma/\nu}$ as a function of field scaling variable $hL^{\beta\delta/\nu}$ for a variety of models. Critical exponents $\gamma$, $\nu$, $\beta$, and $\delta$ are model-dependant. Colored lines and points depict values as measured by the extended algorithm. Solid black lines show a plot of $f(x)=x^{2/\delta}$ for each model. } \label{fig:cluster_scaling} \end{figure*} We have taken several disparate extensions of cluster methods to models in an external field and generalized them to any model of a broad class. This new algorithm has an elegant statement that involves the introduction of not a ghost spin, but a ghost transformation. We provided evidence that extensions deriving from this method are the natural way to extend cluster methods tithe presence of a field, in the sense that it appears to reproduce the scaling of the dynamics in a field that would be expected from renormalization group predictions. In addition to uniting several extensions of cluster methods under a single description, our approach allows the application of fields not possible under prior methods. Instead of simply applying a spin-like field, this method allows for the application of \emph{arbitrary functions} of the spins. For instance, theoretical predictions for the effect of symmetry-breaking perturbations on spin models can be tested numerically \cite{jose_renormalization_1977} \cite{blankschtein_fluctuation-induced_1982,bruce_coupled_1975,manuel_carmona_$n$-component_2000}. \begin{acknowledgments} \end{acknowledgments} \bibliography{monte-carlo} \end{document}