path: root/monte-carlo.tex

%
%  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}