lots of changes

1 files changed, 306 insertions, 433 deletions

diff --git a/monte-carlo.tex b/monte-carlo.tex
index 7133461..929fc65 100644
--- a/monte-carlo.tex
+++ b/monte-carlo.tex
@@ -6,7 +6,7 @@
 \documentclass[aps,prl,reprint]{revtex4-1}
 
 \usepackage[utf8]{inputenc}
-\usepackage{amsmath,amssymb,latexsym,mathtools,xifthen}
+\usepackage{amsmath,amssymb,latexsym,mathtools,xifthen,upgreek}
 \usepackage{algcompatible}
 \usepackage{algorithm}
 
@@ -23,6 +23,7 @@
 \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
@@ -30,6 +31,8 @@
 \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}
@@ -74,6 +77,9 @@
 \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
@@ -82,7 +88,7 @@
 
 \begin{document}
 
-\title{An efficient Wolff algorithm in an external magnetic field}
+\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}
@@ -90,472 +96,339 @@
 \date\today
 
 \begin{abstract}
-  We introduce an extension of the Wolff algorithm that preforms efficiently
-  in an external magnetic field. Near the Ising critical point, the
-  correlation time of our algorithm has a conventional scaling form that
-  reduces to that of the Wolff algorithm at zero field and becomes more efficient
-  at any nonzero field. As an application, we
-  directly measure scaling functions of observables in the metastable state of
-  the 2D Ising model.
 \end{abstract}
 
 \maketitle
 
-The Ising model is a simple model of a magnet comprised of discrete locally
-interacting one-component spins. Like most systems in statistical mechanics,
-calculation of its ensemble properties cannot be done explicitly and is often
-performed using Monte Carlo methods. Near its continuous phase transition,
-divergent correlation length leads to divergent correlation time in any
-locally-updating algorithm, hampering computation.  With no external field,
-this critical slowing-down is largely alleviated by cluster algorithms---the
-most efficient of which is the Wolff algorithm---whose dynamics are nonlocal
-since each step flips a cluster of spins whose average size diverges with the
-correlation length. While less efficient cluster algorithms, like
-Swendsen--Wang, have been modified to perform in nonzero field, Wolff
-only works at zero field. We describe an extension of the Wolff
-algorithm that works in arbitrary external field while preserving Wolff's efficiency throughout the entire temperature--field parameter
-space.
-
-The Wolff algorithm works by first choosing a random spin and adding it to an
-empty cluster. Every neighbor of that spin that is pointed in the same
-direction as the spin is also added to the cluster with probability
-$1-e^{-2\beta J}$, where $\beta=1/T$ is inverse temperature and $J$ is the
-coupling between sites. This process is repeated for the neighbors of every
-spin added to the cluster. When all sites surrounding the cluster have been
-exhausted, the cluster is flipped. Our algorithm is a simple extension of
-this. An extra spin---often referred to as the ``ghost spin''---is introduced
-and made a nearest neighbor of all others with coupling $|H|$, the magnitude
-of the external field.  The traditional Wolff algorithm is then preformed on
-this new extended lattice exactly as described above, with the extra spin
-treated no differently from any others, i.e., allowed to be added to clusters
-and subsequently flipped. Observables in the original system can be exactly
-estimated using the new one by a simple correspondence.
-
-This paper is divided into three sections. First, the Ising model (and our
-notation for it) is introduced, along with extant Monte Carlo algorithms.
-Second, we introduce our algorithm in detail and compare its efficiency with
-the existing ones. Finally, we use this new algorithm to directly measure
-critical scaling functions for observables of the 2D Ising model in its
-metastable state.
-
-\section{Introduction}
-
-Consider an undirected graph $G=(V,E)$ describing a network of interacting spins. The set
-of vertices $V=\{1,\ldots,n\}$ enumerates the sites of the network, and the
-set of edges $E$ describes connections between neighboring sites. On each site
-is a spin that can take any value from the set $S=\{-1,1\}$. The state of the
-system is described by a function $s:V\to S$, leading to a configuration space
-of all possible states $S^n=S\times\cdots\times S$. The Hamiltonian
-$\H:S^n\to\R$ gives the energy of a particular state $s\in S^n$ and is defined
-by
+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 and $X$ an $R$-set, 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 $O(n)$ on $S^{n-1}$ in the $n$-component
+model. 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{loos1969symmetric} 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(s)=-\sum_{\{i,j\}\in E}J_{ij}s_is_j-HM(s),
-  \label{eq:ham}
+  \H(\vec s)=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!Z(s_i,s_j)-\sum_{i\in V}B(s_i),
 \]
-where $H$ is the external magnetic field, $J:E\to\R$ gives the coupling
-between spins on neighboring sites and
-$M:S^n\to\R$ is the magnetization of the system defined for a state
-$s\in S^n$ by
+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 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}
+          & Spin states ($X$) & Symmetry ($R$)  &  Action ($g\cdot s$)  &
+    Coupling ($Z(s,t)$) & Field ($B(s)$) \\
+    \hline\hline
+    Ising         & $\{-1,1\}$    & $\Z/2\Z$                 &  $0\cdot
+    s\mapsto s$, $1\cdot s\mapsto -s$ & $st$ & $Hs$ \\
+    $n$-component & $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)$\\
+    Roughening & $\mathbb Z$ & $D_\inf$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$
+    & $(s-t)^2$ & $Hs^2$\\
+  \end{tabular}
+  \caption{Some models.}
+  \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
 \[
-  M(s)=\sum_{i\in V}s_i.
+\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)}
 \]
-For the purpose of this study, we will only be considering ferromagnetic
-systems where the function $J$ is nonnegative.
-
-An observable of the system is a function $A:S^n\to\R$ depending on the
-system's state. Both the Hamiltonian and magnetization defined above are
-observables. Assuming the ergodic hypothesis holds, the expected value $\avg
-A$ of
-any observable $A$ is defined by its average over every state $s$ in the
-configuration space $S^n$ weighted by the Boltzmann factor $e^{-\beta\H(s)}$,
-or
+
+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 $C\subseteq E$ is related to the probability of the reverse process $P(\vec
+s'\to\vec s)$ by
+\begin{widetext}
 \[
-  \avg A=\frac1Z\sum_{s\in S^n}e^{-\beta\H(s)}A(s),
-  \label{eq:avg.obs}
+  \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}
 \]
-where the partition function $Z$ is defined by
+\end{widetext}
+which satisfies detailed balance so long as $p_r(s,t)=p_{r^{-1}}(s,t)$.
+Therefore, for the algorithm to function one must choose rotations at random
+only from those group elements in $R$ that act by idempotents. Usually, as is
+the case in all the examples here, the action is natural enough that this is
+the same as only choosing group elements of order two to act on your spins.
+Ergodicity is achieved if the subset of symmetry rotations that obey this
+property allow any spin state to access any other under composition. This is
+true for all the
+models in Table \ref{table:models}.
+
+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\}$
+and
 \[
-  Z=\sum_{s\in S^n}e^{-\beta\H(s)}
+  \tilde E=E\cup\big\{\{0,i\}\mid i\in V\big\}.
 \]
-and gives the correct normalization for the weighted sum.
-Unfortunately, the sum over configurations in \eqref{eq:avg.obs} are intractable for all but
-for very small systems. Therefore expectation values are usually
-approximated by some means. Monte Carlo methods are a common way of
-accomplishing this. These methods sample states $s$ from the configuration space
-$S^n$ according to the Boltzmann distribution $e^{-\beta\H(s)}$ so that
-averages of observables made using their incomplete samples asymptotically approach the
-true expected value.
-
-The Metropolis--Hastings algorithm
-\cite{metropolis1953equation,hastings1970monte} is very popular for systems in statistical
-physics. A starting state $s\in S^n$ is randomly perturbed to the state $s'$,
-usually by flipping one spin. The change in energy $\Delta\H=\H(s')-\H(s)$ due
-to the perturbation is then computed. If the change is negative the perturbed
-state $s'$ is accepted as the new state. Otherwise the perturbed state
-is accepted with probability $e^{-\beta\Delta\H}$. This process is repeated
-indefinitely and a sample of states is made by sampling the state $s$
-between iterations sufficiently separated that the successively sampled states
-are uncorrelated, e.g., at separations larger than the correlation time $\tau$.  Metropolis--Hastings is very general, but unless the
-perturbations are very carefully chosen the algorithm suffers in regimes where
-large correlations are present in the system, for instance near continuous
-phase transitions. Here the algorithm suffers from what is known as critical
-slowing-down, where likely states consist of large correlated clusters that
-take many perturbations to move between in configuration space.
-
-%\begin{figure}
-%  \begin{algorithm}[H]
-%    \begin{algorithmic}
-%      \REQUIRE $s\in S^n$
-%      \STATE $s'\gets$ \texttt{Perturb}($s$)
-%      \STATE $\Delta\H\gets\H(s')-\H(s)$
-%      \IF {$\exp(-\beta\Delta\H)>$ \texttt{UniformRandom}$(0,1)$}
-%        \STATE $s\gets s'$
-%      \ENDIF
-%    \end{algorithmic}
-%    \caption{Metropolis--Hastings}
-%    \label{alg:met-hast}
-%  \end{algorithm}
-%\end{figure}
-%
-The Wolff algorithm \cite{wolff1989collective} solves many of these problems, but only at zero external
-field, $H=0$. This algorithm solves the problem of critical slowing-down by
-flipping carefully-constructed clusters of spins at once in a way that samples
-high-correlated states quickly while also always accepting prospective states.
-A random site $i$ is selected from the graph and its spin is flipped. Each of the
-site's neighbors $j$ is also flipped with probability $1-e^{-2\beta J_{ij}}$
-if doing so would lower the energy of the bond $i$--$j$. The process is
-repeated with every neighbor that was flipped. While this algorithm
-successfully addresses the problems of critical slowing down at zero field, it
-doesn't directly apply at nonzero field. A common extension is the hybrid
-Wolff--Metropolis, where a Wolff cluster is formed using the ordinary formula,
-then accepted or rejected Metropolis-style based on the change in energy
-\emph{relative to the external field}, $H\Delta M$. Because each spin has the
-same coupling to the external field, a strictly more efficient
-but exactly equivalent version of this hybrid is made by distributing the
-cluster acceptance or rejection across the spins in the cluster as they are
-added. In this version, the cluster is abandoned with probability $1-e^{-2\beta
-H}$ every time a spin is added to it.
-
-$z=0.29(1)$ \cite{wolff1989comparison,liu2014dynamic} $z=0.35(1)$ for
-Swendsen--Wang \cite{swendsen1987nonuniversal}
-
-%\begin{figure}
-%  \begin{algorithm}[H]
-%    \begin{algorithmic}
-%      \REQUIRE $s\in\G^n$
-%      \STATE \textbf{let} $q$ be an empty stack
-%      \STATE \textbf{let} $i_0\in V$
-%      \STATE $q.\mathtt{push}(i_0)$
-%      \STATE $\sigma\gets s_{i_0}$
-%      \WHILE {$q$ is not empty}
-%        \STATE $i\gets q.\mathtt{pop}$
-%        \IF {$s_i=\sigma$}
-%          \FORALL {$j$ such that $\{i,j\}\in E$}
-%            \IF {$1-\exp(-2\beta
-%              J_{ij}s_is_j)>\mathop{\mathtt{UniformRandom}}(0,1)$}
-%          \STATE $q.\mathtt{push}(j)$
-%            \ENDIF
-%          \ENDFOR
-%          \STATE $s_i\gets-s_i$
-%        \ENDIF
-%      \ENDWHILE
-%    \end{algorithmic}
-%    \caption{Wolff (Zero-Field)}
-%    \label{alg:wolff}
-%  \end{algorithm}
-%\end{figure}
-%
-%
-%
-%\begin{figure}
-%  \begin{algorithm}[H]
-%    \begin{algorithmic}
-%      \REQUIRE $s\in\G^n$
-%      \STATE $s'\gets s$
-%      \STATE \textbf{let} $q$ be an empty stack
-%      \STATE \textbf{let} $i_0\in V$
-%      \STATE $q.\mathtt{push}(i_0)$
-%      \STATE $\sigma\gets s_{i_0}'$
-%      \STATE \texttt{completed} $\gets$ \textbf{true}
-%      \WHILE {$q$ is not empty}
-%        \STATE $i\gets q.\mathtt{pop}$
-%        \IF {$s_i'=\sigma$}
-%          \FORALL {$j$ such that $\{i,j\}\in E$}
-%          \IF {$1-\exp(-2\beta
-%            J_{ij}s_i's_j')>\mathop{\mathtt{UniformRandom}}(0,1)$}
-%          \STATE $q.\mathtt{push}(j)$
-%            \ENDIF
-%          \ENDFOR
-%          \STATE $s_i'\gets-s_i'$
-%          \STATE $q.\mathtt{push}(i)$
-%          \IF {$1-\exp(-2\beta\sigma H)>\mathop{\mathtt{UniformRandom}}(0,1)$}
-%          \STATE \texttt{completed} $\gets$ \textbf{false}
-%            \STATE \textbf{break}
-%          \ENDIF
-%        \ENDIF
-%      \ENDWHILE
-%      \IF {completed}
-%        $s\gets s'$
-%      \ENDIF
-%    \end{algorithmic}
-%    \caption{Hybrid Wolff/Metropolis--Hastings}
-%    \label{alg:h.wolff}
-%  \end{algorithm}
-%\end{figure}
-
-\section{Cluster-flip in a Field}
-
-Consider the new graph $\tilde G=(\tilde V,\tilde E)$ defined from $G$ by
-\begin{align}
-  \tilde V&=0\cup V\\
-  \tilde E&=E\cup\{\{0,i\}\mid i\in V\},
-\end{align}
-or by adding a zeroth vertex and edges from every other vertex to the new
-one. Thee spin at site zero is often known as a ``ghost spin.'' The network of spins
-described by this graph now have states that occupy a configuration space
-$S^{n+1}$. Extend the coupling between spins $\tilde J:\tilde E\to\R$ by
+We have introduced a new site to the lattice that neighbors every other site.
+Instead of assigning this site a spin whose value comes from the set $X$, we
+will assign it values $s_0\in R$, symmetry group elements, 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
 \[
-  \tilde J(e)=\begin{cases}
-    J(e) & e\in E\\
-    |H| & \text{otherwise},
-  \end{cases}
+\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}
 \]
-so that each spin is coupled to every other spin the same way they were
-before, and coupled to the new spin with strength $|H|$, the magnitude of the
-external field. Now define a
-Hamiltonian function 
-$\tilde\H:S^{n+1}\to\R$ on this new, larger configuration space defined by
+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\H(s)=-\sum_{\{i,j\}\in\tilde E}\tilde J_{ij}s_is_j.
-  \label{eq:new.h.simple}
+  \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}
 \]
-This new Hamiltonian resembles the old one, but is comprised only of
-spin--spin interactions with no external field. However, by changing the terms
-considered in the sum we may equivalently write
+Note that this modified coupling is invariant under the action of group
+elements: for any $r,s_0\in R$ and $s\in X$,
 \[
-  \tilde\H(s)=-\sum_{\{i,j\}\in E}J_{ij}s_is_j-H\tilde M(s)
+\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}
 \]
-where the new magnetization $\tilde M:S^{n+1}\to\R$ is defined for $s\in
-S^{n+1}$ by
+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
 \[
-  \begin{aligned}
-    \tilde M(s)
-    &=\sgn(H)s_0\sum_{i\in V}s_i\\
-    &=M(s_0\sgn(H)(s_1,\ldots,s_n)).
-  \end{aligned}
+  \tilde A(s_0,\vec s)=A(s_0^{-1}\cdot\vec s)
 \]
-In fact, any observable $A$ of the original system can be written as an
-observable $\tilde A$ of the new system by defining
+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
 \[
-  \tilde A(s)=A(s_0\sgn(H)(s_1,\ldots,s_n))
+\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}
 \]
-such that $\avg{\tilde A}=\avg A$. This can be seen readily by using the
-symmetry $\tilde\H(-s)=\tilde\H(s)$ of the Hamiltonian
-\eqref{eq:new.h.simple}, the invariance of configuration space sums under
-negation of their summand, and the fact that $\tilde\H(1,s)=\H(s)$ for any $s\in
-S^n$. Notice in particular that this is true for the Hamiltonian $\tilde\H$ as
-well.
-
-Our new spin system with an additional field is simply a ferromagnetic spin system in the absence of an external field.
-Therefore, the Wolff algorithm can be applied to it with absolutely no
-modifications. Since there is an exact correspondence between the expectation
-values of our ordinary spin system in a field and their appropriately defined
-counterparts in our new system, sampling Monte Carlo simulations of our new
-system allows us to estimate the expectation values for the old. This ``new''
-algorithm, if you can call it that, is shown in Algorithm
-\ref{alg:wolff-field}.
+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.
+
+\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.
+
+\subsection{The $n$-component Model}
+
+In the $n$-component 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.
+
+\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{Roughening}
+
+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|$, though it may also be any function of the absolute difference.
+Because randomly 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.
+
 
 %\begin{figure}
-%  \begin{algorithm}[H]
-%    \begin{algorithmic}
-%      \REQUIRE $s\in S^n$
-%      \STATE \textbf{let} $q$ be an empty stack
-%      \STATE \textbf{let} $i_0\in V$
-%      \STATE $q.\mathtt{push}(i_0)$
-%      \STATE $\sigma\gets s_{i_0}$
-%      \WHILE {$q$ is not empty}
-%        \STATE $i\gets q.\mathtt{pop}$
-%        \IF {$s_i=\sigma$}
-%          \FORALL {$j$ such that $\{i,j\}\in\tilde E$}
-%            \IF {$1-\exp(-2\beta \tilde
-%              J_{ij}s_is_j)>\mathop{\mathtt{UniformRandom}}(0,1)$}
-%          \STATE $q.\mathtt{push}(j)$
-%            \ENDIF
-%          \ENDFOR
-%          \STATE $s_i\gets-s_i$
-%        \ENDIF
-%      \ENDWHILE
-%    \end{algorithmic}
-%    \caption{Wolff (Nonzero Field)}
-%    \label{alg:wolff-field}
-%  \end{algorithm}
+%  \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}
-%
-
-At sufficiently small $|H|$ both our algorithm and the hybrid Wolff--Metropolis reduce
-to the ordinary Wolff algorithm, and have its runtime properties. At very
-large $|H|$, the hybrid
-Wolff--Metropolis behaves exactly like Metropolis, where almost only one spin
-is ever flipped at a time with probability $\sim e^{-2\beta H}$, since the
-energy is dominated by contributions from the field.
-
-We measured the autocorrelation time $\tau$ of the internal energy $\H$ of a
-square-lattice Ising model ($J=1$ and $E$ is the set of nearest neighbor
-pairs)
-resulting from using each of these three algorithms at various
-fields and temperatures. This was done using an initial convex sequence
-estimator, as described in \cite{geyer1992practical}. Time was
-measured as ``the number of spins that the algorithm has attempted to flip.''
-For example, every Metropolis--Hastings step takes unit time, every Wolff step takes
-time equal to the number of spins in the flipping cluster, and every hybrid
-Wolff/Metropolis--Hastings step takes time equal to the number of spins in the
-cluster the step attempted to flip, whether it succeeds or not. The
-autocorrelation time as a function of reduced field $\beta H$ is shown in
-Fig.~\ref{fig:correlation}. By this measure of time, the modified Wolff
-algorithm thermalizes faster than both the hybrid Wolff/Metropolis and
-ordinary Metropolis--Hastings in all regimes. At high field ($\beta H\sim1$)
-the new algorithm has an autocorrelation time of roughly $n$, while both
-Metropolis--Hastings and the hybrid algorithm have an autocorrelation time of
-$2n$. At low temperature, neither the new algorithm nor Metropolis-Hastings
-have an autocorrelation time that depends on field, but the hybrid method
-performs very poorly at intermediate values of field in a way that becomes
-worse with system size, scaling worse than the hybrid algorithm does even it
-the critical temperature. At high temperature, all three algorithms perform
-similarly and in a way that is independent of system size. At the critical
-temperature, both the Metropolis--Hastings and hybrid algorithms experience
-critical slowing down in a similar way, whereas the new algorithm slows down
-in a far reduced way, always performing more efficiently at nonzero field than
-the zero-field Wolff algorithm.
 
 \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}
+  \input{fig_correlation_collapse-hL}
+  \caption{The correlation time $\tau$ as a function of the renormalization
+    invarient $hL^{-\beta\delta/\nu}$ for the $N=L\times L$ square lattice
+    Ising model for $L=8$, $16$, $32$, $64$, $128$, and $256$. $z=0.3$
+  }
 \end{figure}
 
 \begin{figure}
   \centering
-  \input{fig_correlation_collapse-hL}
+  \input{fig_correlation-temp}
+  \caption{The correlation time $\tau$ as a function of the renormalization
+    invarient $ht^{-\beta\delta}$ for the $N=128\times 128$ square lattice
+    Ising model. $z=0.3$
+  }
 \end{figure}
 
-Our algorithm for the Ising model in a field can be
-generalized to run on the $q$-spin Potts or $O(n)$ models in exactly the same
-way as the conventional Wolff algorithm.
-
-\section{Magnetization Estimator}
-
-At any size, the ensemble average magnetization $\avg M$ is identically zero
-at zero field. However, this is not what is observed at low temperature. This
-is because in the low temperature phase the ergodic hypothesis---that the
-time-average value of observables is equal to their ensemble average---is
-violated. As the system size grows the likelihood of a fluctuation in any
-reasonable dynamics that flips the magnetization from one direction to the
-other becomes vanishingly small, and therefore it is inappropriate to estimate
-expected values in the low temperature phase by averaging over the whole
-configuration space. Instead, values must be estimated by averaging over the
-portion of configuration space that is accessible to the dynamics.
-
-For finite size systems, like any we would simulate, dynamics at zero field or
-even small nonzero field do allow the whole configuration space to be
-explored. However, people usually want to use the results from finite size
-systems to estimate the expected values in the thermodynamic limit, where this
-is no longer true. At zero field, for instance, it is common practice to use
-$\avg{|M|}$ to estimate the expected value for the magnetization instead of
-$\avg M$. But what to do at finite field? Is this approach justified?
-
-The provenance of the restricted configuration space in the thermodynamic
-limit is the fact that, at criticality or for nonzero $H$ under the critical
-temperature, an infinite cluster forms and can never be flipped by the
-dynamics. In a finite system there is no infinite cluster, and even the
-largest cluster has a nonzero probability of being flipped. However, we still
-want to be able to use finite systems to estimate the quantities of infinite
-ones. We approach this by thinking of a finite system as a piece of an
-infinite one. Clearly, time-averages of quantities in the chunk are equal to
-time-averages of the whole system, simply with much slower convergence.
-However, simulation of finite systems present a problem: in each sample, we
-are given a snapshot of a chunk of the system but do not know which direction
-the infinite cluster in the surrounding infinite system is pointing.
-Therefore, we must make a guess such that, in the limit of larger and larger
-system size, the probability that our guess is wrong and that the infinite
-cluster is pointed in a direction opposite the one we expected approaches
-zero. To estimate the equilibrium values of quantities in the infinite system,
-we only admit values resulting from states whose corresponding infinite
-cluster is expected to point in the direction of our external field. If the
-external field is zero, we choose one direction at random, say positive.
-
-Snapshots of the chunk what we suspect are likely to be a piece of an infinite
-cluster facing opposite the external field are discarded for the purpose of
-computing equilibrium values of the system, but still represent something
-physical. When the model's state contains an infinite cluster oriented against
-the external field, it is known as \emph{metastable}. Therefore in a finite
-simulation one can make convergent estimates of quantities in both the
-equilibrium and metastable states of the infinite system by dividing sampled
-states into two sets: those that are likely snapshots of a finite block in an
-infinite system whose infinite cluster is oriented with the field, and those
-whose infinite cluster is oriented against the field. All we must do now is
-specify how to distinguish between likely oriented with and likely oriented
-against. Others have used the direction of the largest cluster in the finite
-system to make this distinction. Far from the critical point in the
-low-temperature phase this does a clear and consistent job. However, we would
-like to study systems very near the critical point. We will take the likely direction of the
-infinite cluster to be given by \emph{the direction in which the system is
-magnetized}. Clearly, as our finite chunk grows larger and larger, the
-probability that its magnetization is different from the direction of the
-infinite cluster goes to zero. In the thermodynamic limite
-
-Since, in the thermodynamic limit expected values are given by an average
-over a restricted configuration space, we can estimate those expected values
-at finite size by making the same restriction. Defining the reduced
-configuration spaces
-\begin{align}
-  S_\e^n&=\{s\in S^n\mid \sgn(H)M(s)>0\}
-  \\
-  S_\m^n&=\{s\in S^n\mid \sgn(H)M(s)<0\}
-  \\
-  S_0^n&=\{s\in S^n\mid \sgn(H)M(s)=0\}
-\end{align}
-where
-\[
-  \sgn(H)=\begin{cases}1&H\geq0\\-1&H<0.\end{cases}
-\]
-Clearly, $S^n=S^n_\e\cup S^n_\m\cup S^n_0$, which implies that
-$Z=Z_\e+Z_\m+Z_0$. Since $S^n_\e=\{-s\mid s\in S^n_\m\}$, $|S^n_\e|=|S^n_\m|$.
-At zero field, $\H(s)=\H(-s)$, and therefore $Z_\e=Z_\m$.  This can be used to
-show the near-equivalence (at zero field) of taking expectation values
-$\avg{|M|}$ of the absolute value of the magnetization and taking expectation
-values $\avg M_\e$ of the magnetization on a reduced configuration space,
-since
-\begin{align}
-  \avg{|M|}=\frac1{1+Z_0/2Z_\e}\eavg M
-\end{align}
-At infinite temperature, $Z_0/Z_\e\simeq n^{-1/2}\sim L^{-1}$ for large $L$,
-$n$. At any finite temperature, especially in the ferromagnetic phase,
-$Z_0\ll Z_\e$ in a much more extreme way.
-
-If the ensemble average over only positive magnetizations can be said to
-converge to the equilibrium state of the Ising model in the thermodynamic
-limit, what of the average over only \emph{negative} magnetizations, or the
-space $S^n_\m$,
-ordinarily unphysical?  This is, in one sense, precisely the definition of the
-metastable state of the Ising model. Expected values of observables in the
-metastable state can therefore be computed by considering this reduced
-ensemble.
-
-How, in practice, are these samplings over reduced ensembles preformed? There
-may be efficient algorithms for doing Monte Carlo sampling inside a particular
-reduced ensemble, but for finite-size systems with any of the algorithms
-described here the whole configuration space is available. Any of these can be
-used, however, to simultaneously sample all reduced configuration spaces by
-allowing them to sample the whole space and then only add a contribution to
-a given average if the associated state is in the reduced space of interest.
-
 \begin{acknowledgments}
-  Thanks!
 \end{acknowledgments}
 
 \bibliography{monte-carlo}