%
%  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{Accelerating Monte Carlo: Wolff in arbitrary external fields}
\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 natural way to extend celebrated spin-cluster Monte Carlo
  algorithms for fast thermal lattice simulations at criticality, like Wolff, to
  systems in arbitrary fields. The method relies on the generalization of the
  `ghost spin' representation to one with a `ghost transformation' that
  restores invariance to spin symmetries at the cost of an extra degree of
  freedom. The ordinary cluster-building process can then be run on the new
  representation. For several canonical systems, we show that this extension
  preserves the scaling of accelerated dynamics in the absence of a field.
\end{abstract}

\maketitle

Lattice models are important in the study of statistical physics and phase
transitions. Rarely exactly solvable, they are typically studied by
approximate and numerical methods. Monte Carlo techniques are a common way of
doing this, approximating thermodynamic quantities by sampling the
distribution of systems states. These Monte Carlo algorithms are better the
faster they arrive at a statistically independent sample. This typically
becomes a problem near critical points, where critical slowing down
\cite{wolff_critical_1990} results in power-law divergences of dynamic
timescales. Celebrated cluster algorithms largely addressed this in the absence of symmetry-breaking 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}.
These result in relatively small dynamic exponents for many spin
systems \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
symmetry operations on the spins. Some success has been made in extending these
algorithms to systems in certain external fields by adding a `ghost site'
\cite{coniglio_exact_1989} that returns
global rotation invariance to spin Hamiltonians at the cost of an extra degree
of freedom, allowing the method to be used in a subcategory of interesting
fields \cite{alexandrowicz_swendsen-wang_1989, wang_clusters_1989, ray_metastability_1990}. Static
fields have also been applied by including a separate metropolis or heat bath
update step after cluster formation \cite{destri_swendsen-wang_1992,
lauwers_critical_1989}, and other categories of fields have
been applied using replica methods
\cite{redner_graphical_1998,chayes_graphical_1998,machta_replica-exchange_2000}. We show that the scaling of
correlation time near the critical point of several models suggests that the
`ghost' 
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 spin systems, \emph{arbitrary} external fields can be treated
using cluster methods. Rather than the introduction of a `ghost spin,\!' our
representation relies on introducing a `ghost transformation.\!'

We will pose the problem in a general way, but several specific examples can
be found in Table~\ref{table:models} for concreteness. 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 their set of isometries.  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 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. The formal results of this
paper (that the algorithm obeys detailed balance and ergodicity) hold equally
well for these cases, but we will drop the additional index notation for
clarity. Statements about efficiency may not.

\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 the $\mathrm O(2)$ model \cite{jose_renormalization_1977}.}
  \label{table:models}
\end{table*}

The goal of statistical mechanics 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 Boltzmann probability 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 using 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 of several related cluster algorithms can be described for this
system, we will focus on the Wolff algorithm \cite{wolff_collective_1989}. In
the absence of an external field, e.g., B(s)=0, 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. 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
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)}\\
  &\quad=\prod_{\{i,j\}\in\partial
  C}e^{\beta(Z(r\cdot s_i,s_j)-Z(s_i,s_j))}
  =\frac{p_r(s_i,s_j)}{p_{r}(s_i,s_j)}\frac{e^{-\beta\H(\vec
  s)}}{e^{-\beta\H(\vec s')}},
\end{aligned}
\]
%\end{widetext}
whence detailed balance is also satisfied. 

This algorithm relies on the fact that the coupling $Z$ depends only on
relative orientation of the spins---global reorientations do not affect the
Hamiltonian. The external field $B$ breaks this symmetry. However, it can be
restored. Define a new graph $\tilde G=(\tilde V,\tilde E)$, where $\tilde
V=\{0,1,\ldots,N\}$ adds the 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 the ghost site a spin whose value
comes from $X$, we assign it values in the symmetry group $s_0\in R$, so that
the configuration space of the new model is $R\times X^N$. We introduce the
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}
\]
The 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}\cdot s)
  =\tilde Z(s_0,s)
\end{aligned}
\]
The invariance of $\tilde Z$ to rotations given other arguments follows from
the invariance properties of $Z$.

We have produced a system incorporating the field function $B$ whose
Hamiltonian is invariant under global rotations, but how does it relate to our
old system, whose properties we actually want to measure? If $A:X^N\to\R$ is
an observable of the original system, we 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 invariant under global rotations, it follows
that for any $g\in R$, $\tilde\H(g,g\cdot\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}
\]
Using this equivalence, 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 whose value 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 might be considered the `ghost
transformation' representation.

Several specific examples from Table~\ref{table:models} are described in the
following.

\emph{The Ising model.} In the Ising model spins are drawn from the set $\{1,-1\}$. Its symmetry group
is $C_2$, the cyclic group on two elements, which can be conveniently
represented by a multiplicative group with elements $\{1,-1\}$, exactly the
same as the spins themselves. The only nontrivial element is of order two.
Since the symmetry group and the spins are described by the same elements,
performing the algorithm on the Ising model in a field is fully described by
just using the `ghost spin' representation.  This algorithm or algorithms
based on the same decomposition of the Hamiltonian have been applied
by several researchers \cite{alexandrowicz_swendsen-wang_1989,
wang_clusters_1989, ray_metastability_1990}. The algorithm has been
implemented by one of the authors in an existing interactive Ising
simulator at \texttt{https://mattbierbaum.github.io/ising.js} \cite{bierbaum_ising.js_nodate}.

\emph{The $\mathrm O(n)$ model.} In the $\mathrm O(n)$ model spins are described by vectors on the
$(n-1)$-sphere $S^{n-1}$. Its symmetry group is $O(n)$, $n\times n$ orthogonal
matrices, which act on the spins by matrix multiplication. The elements of
$O(n)$ of order two are reflections about hyperplanes through the origin and
$\pi$ rotations about any axis through the origin. Since the former generate
the entire group, reflections alone suffice to provide ergodicity. The `ghost
spin' version of the algorithm has been used to apply a simple vector field to
the $\mathrm O(3)$ model \cite{dimitrovic_finite-size_1991}. Other fields of
interest include $(n+1)$-dimensional spherical harmonics
\cite{jose_renormalization_1977} and cubic fields
\cite{bruce_coupled_1975,blankschtein_fluctuation-induced_1982}, which can be
applied with the new method. The method is
quickly generalized to spins whose symmetry groups other compact Lie groups

\emph{The Potts \& clock models.} In both the $q$-state Potts and clock models spins are described by elements
of $\Z/q\Z$, the set of integers modulo $q$. Its symmetry group 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 by
$2\pi n/q$, and the element $s_n$ represents a reflection composed with the
rotation $r_n$. The group acts on spins by permutation: $r_n\cdot m={n+m}\pmod
q$ and $s_n\cdot m={-(n+m)}\pmod q$. This is 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$ of order 2 are all reflections and $r_{q/2}$ if
$q$ is even, though the former can generate the latter. While reflections do
not necessarily generate the entire group, their action on $\Z/q\Z$ is
transitive and therefore the algorithm is ergodic.

\emph{Roughening models.} Though not often thought of as a spin model, roughening of surfaces can be
described in this framework. Spins are described by integers $\Z$ and their
symmetry group is the infinite dihedral group $D_\infty=\{r_i,s_i\mid
i\in\Z\}$, whose action on the spin $j\in\Z$ is given by $r_i\cdot j=i+j$ and
$s_i\cdot j=-i-j$. The elements of order two are reflections $s_i$, whose
action on $\Z$ is transitive. The coupling can be any function of the absolute
difference $|i-j|$.  Because random choice of reflection will almost always
result in energy changes so large that the whole system is flipped, it is
better to select random reflections about integers close to the average state
of the system.  A variant of the algorithm has been applied without a field
\cite{evertz_stochastic_1991}.


No algorithm is worthwhile if it doesn't run efficiently. This algorithm,
being an extension of the Wolff algorithm into a new domain, should be
considered successful if it likewise extends the efficiency of the Wolff
algorithm into that domain. Some systems are not efficient under Wolff, and we
don't expect this extension to help them. For instance, Ising models with
random fields or bonds technically can be treated with Wolff
\cite{dotsenko_cluster_1991}, but it is not efficient because the clusters
formed do scale naturally with the correlation length
\cite{rieger_monte_1995,redner_graphical_1998}. Other approaches, like replica methods, should
be relied on instead
\cite{redner_graphical_1998,chayes_graphical_1998,machta_replica-exchange_2000}.

At a critical point, correlation time $\tau$ scales with system size
$L=N^{-D}$ as $\tau\sim L^z$. Cluster algorithms are celebrated for their
small dynamic exponents $z$. In the vicinity of an ordinary critical point,
the renormalization group predicts scaling behavior for the correlation time
as a function of temperature $t$ and field $h$ of the form
\[
  \tau=h^{-z\nu/\beta\delta}\mathcal T(ht^{-\beta\delta},hL^{\beta\delta/\nu}).
\]
If a given dynamics for a system at zero field results in scaling like $L^z$,
one should expect its natural extension in the presence of a field to scale
roughly like $h^{-z\nu/\beta\delta}$ and collapse appropriately as a function
of $hL^{\beta\delta/\nu}$.  We measured the autocorrelation time for the $D=2$
square-lattice model at a variety of system sizes, temperatures, and fields
$B(s)=hs/\beta$ using standard methods \cite{geyer_practical_1992}. The
resulting scaling behavior, plotted in
Fig.~\ref{fig:correlation_time-collapse}, is indeed consistent with an
extension to finite field of the behavior at zero field.

\begin{figure}
  \centering
  \input{fig_correlation_collapse-hL}
  \caption{Collapse of the correlation time $\tau$ of the 2D square lattice
    Ising model 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}$. The exponent
    $z=0.30$ is taken from recent measurements at zero field
  \cite{liu_dynamic_2014}. The solid black line shows a plot of
  $(hL^{\beta\delta/\nu})^{-z\nu/\beta\delta}$.
  }
  \label{fig:correlation_time-collapse}
\end{figure}

Since the formation and flipping of clusters is the hallmark of 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 fact that the 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 Swendsen--Wang 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)\\
    &=L^{\gamma/\nu}g(ht^{-\beta\delta},hL^{\beta\delta/\nu}).
  \end{aligned}
\]

For the Ising model, an additional scaling relation can be written. Since 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},ht^{-\beta\delta})^{-\gamma/\beta\delta}\beta \mathcal
      Y(hL^{\beta\delta/\nu,ht^{-\beta\delta}})\\
      &\hspace{1em}+(hL^{\beta\delta/\nu},ht^{-\beta\delta})^{2/\delta}\mathcal
    M(hL^{\beta\delta/\nu},ht^{-\beta\delta})\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 $g(0,x)\propto x^{2/\delta}$ for each model.
  }
  \label{fig:cluster_scaling}
\end{figure*}

We have taken several disparate extensions of cluster methods to spin models
in an external field and generalized them to work for any model of a broad
class.  The resulting representation involves the introduction of not a ghost
spin, but a ghost transformation. We provided evidence that algorithmic
extensions deriving from this method are the natural way to extend cluster
methods in the presence of a field, in the sense that they appear to reproduce
the scaling of dynamic properties 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, blankschtein_fluctuation-induced_1982,
bruce_coupled_1975, manuel_carmona_$n$-component_2000}.

\begin{acknowledgments}
\end{acknowledgments}

\bibliography{monte-carlo}


\end{document}