added a figure, did a bunch of rewriting

1 files changed, 219 insertions, 175 deletions

diff --git a/monte-carlo.tex b/monte-carlo.tex
index 6509eab..50cd277 100644
--- a/monte-carlo.tex
+++ b/monte-carlo.tex
@@ -7,8 +7,8 @@
 
 \usepackage[utf8]{inputenc}
 \usepackage{amsmath,amssymb,latexsym,mathtools,xifthen}
+\usepackage{algcompatible}
 \usepackage{algorithm}
-\usepackage{algorithmic}
 
 % uncomment to label only equations that are referenced in the text
 %\mathtoolsset{showonlyrefs=true}
@@ -82,8 +82,6 @@
 
 \begin{document}
 
-\algsetup{indent=1em}
-
 \title{A cluster algorithm for the Ising model in an external field}
 \author{Jaron Kent-Dobias}
 \author{James P.~Sethna}
@@ -99,134 +97,120 @@
 
 \section{Introduction}
 
-The Ising model is a simple model of a magnet. It is comprised of $n$ sites
-arranged in a regular lattice. At every site $i$ is a magnetic spin that can
-be oriented up or down, corresponding to a local magnetization $s_i$ with
-value $1$ or $-1$, respectively. Configurations of the spins are uniquely
-described by the local magnetization of every site, 
-so the model's
-configuration space is $\G^n=\G\times\cdots\times\G$ for $\G=\{-1,1\}$
-and states of the model are given by $s=(s_1,\ldots,s_n)\in\G^n$. The
-Ising Hamiltonian $\H:\G^n\to\R$ gives the energy of a particular state
-$s\in\G^n$ by
+The Ising model is a simple model of a magnet comprised of locally interacting
+spins.
+
+Consider an undirected graph $G=(V,E)$ describing a system 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 interacting 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
 \[
-  \H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-HM(s)
+  \H(s)=-\sum_{\{i,j\}\in E}J_{ij}s_is_j-HM(s),
   \label{eq:ham}
 \]
-where $H$ is the external magnetic field, $M:\G^n\to\R$ is the magnetization
-and is defined for $s\in\G^n$ by
-\[
-  M(s)=\sum_{i=1}^ns_i,
-\]
-and the spin--spin coupling $J$ is defined by
-\[
-  J_{ij}=\begin{cases}
-    1 & \text{$i$, $j$ are neighboring sites}\\
-    0 & \text{otherwise.}
-  \end{cases}
-\]
-The quantities of interest in statistical physics are averaged out over
-timescales large compared to the dynamics of the system. If $A:\G^n\to\R$ is a
-physical quantity that depends on the model's configuration and $s(t)$ gives
-the configuration as a function of time $t$, then
+where $H$ is the external magnetic field, $J:E\to\R$ gives the coupling between spins on connected sites and
+$M:S^n\to\R$ is the magnetization of the system defined for a state
+$s\in S^n$ by
 \[
-  \bar A=\lim_{T\to\infty}\frac1T\int_0^TA(s(t))\,\dd t
+  M(s)=\sum_{i\in V}s_i.
 \]
-Typically the average value of physical quantities is not actually computed by
-taking a time average. Instead, the ergodic hypothesis is applied, which
-stipulates that the system's dynamics access all states $s\in S^n$ with relative
-duration $e^{-\beta\H(s)}$, so that the time average can be written as the
-ensemble average, or the average over configuration space,
+For the purpose of this study, we will only be considering ferromagnetic
+systems where the function $J$ is nonnegative. All formal results can be
+generalized to the antiferromagnetic or mixed cases, but algorithmic
+efficiency for instance cannot.
+An observable of this 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
 \[
-  \avg A=\frac{\sum_se^{-\beta\H(s)}A(s)}{\sum_se^{-\beta\H(s)}}
-  =\frac1Z\sum_se^{-\beta\H(s)}A(s),
+  \avg A=\frac1Z\sum_{s\in S^n}e^{-\beta\H(s)}A(s),
   \label{eq:avg.obs}
 \]
-where $\beta=1/kT$ is the inverse temperature in energy units and
+where $\beta$ is the inverse temperature and the partition function $Z$ defined by
 \[
-  Z=\sum_se^{-\beta\H(s)}
+  Z=\sum_{s\in S^n}e^{-\beta\H(s)}
 \]
-is the partition function.
-
-Monte Carlo methods attempt to estimate quantities of the form
-\eqref{eq:avg.obs} without doing the often intractable sum explicitly by
-sampling the set of possible states such that a state $s$ is sampled with
-probably $e^{-\beta\H(s)}$, appropriately weighting each state's contribution
-to the sum. A common way of doing this is the Markov-chain Monte Carlo
-approach, which steps through configuration space using a Markov process whose
-equilibrium distribution is Boltzmann. One of the most popular of these is the
-Metropolis--Hastings algorithm, shown in Algorithm \ref{alg:met-hast}. One
-spin in the current state $s$ is randomly chosen and perturbed. The new state
-is accepted if the resulting change in energy is negative, or with probability
-$e^{-\beta\Delta\H}$ if the change is positive. Thermodynamic averages of
-functions of the state can be made by periodically sampling from this loop for
-times much longer than the autocorrelation time of the Markov process.
+gives the correct normalization for the weighted sum.
+
+The sum over configurations in \eqref{eq:avg.obs} are intractable for all but
+for very small systems. Therefore expected values of observables are usually
+approximated by some means. Monte Carlo methods are a common way of
+accomplishing this. These methods sample states $s$ in the configuration space
+$S^n$ according to the Boltzmann distribution $e^{-\beta\H(s)}$, so that
+averages of observables made using their samples asymptotically approach the
+true expected value.
+
+The Metropolis--Hastings algorithm 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. This algorithm is shown schematically in Algorithm
+\ref{alg:met-hast}. 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\G^n$
-      \LOOP
-      \STATE $s'\gets s$
-      \STATE $i \gets\text{random integer between $1$ and $n$}$
-      \STATE $s'_i\gets-s_i$
-      \IF{$e^{-\beta(\H(s')-\H(s))}\geq\text{random number between $0$ and $1$}$}
-      \STATE $s\gets s'$
+      \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
-      \ENDLOOP
     \end{algorithmic}
     \caption{Metropolis--Hastings}
     \label{alg:met-hast}
   \end{algorithm}
 \end{figure}
 
-This algorithm is very generic and works well for many models. However, near
-the critical point of the Ising model it suffers from critical slowing down,
-where the autocorrelation time of the algorithm diverges with proximity to the
-critical point.
-The Wolff algorithm is such a method that is very efficient when
-$H=0$, i.e., in the absence of an external field. It is shown in Algorithm
-\ref{alg:wolff}. Briefly, whole clusters of spins are flipped in a way that
-correctly samples states from the Boltzmann distribution while never rejecting
-a flip attempt. A random spin is selected, then flipped. Each of the spin's
-neighbors that was initially oriented the same direction as the spin is then
-flipped with probability $1-e^{-2\beta}$, and the same is true recursively for
-each of those flipped spins' neighbors.
-
-The Wolff algorithm solves many of the problems of the Metropolis--Hastings
-algorithm, but only in the absence of a field, since it depends on the natural
-symmetry $\H(s)=\H(-s)$ which is true only when $H=0$. The Wolff algorithm
-can be extended to work in a field by combining it with aspects of the
-Metropolis--Hastings algorithm in what we will refer to as the hybrid
-Wolff/Metropolis--Hastings. In each step of this algorithm, a Wolff-style
-cluster is constructed exactly as above. Instead of moving on, however, the
-system is reverted to its previous state with probability
-$\max\{0,1-e^{\beta\Delta M H}\}$, where $\Delta M$ is the change in
-magnetization due to the cluster flip. This is precisely a
-Metropolis--Hastings step where the perturbation of $s$ to $s'$ is the Wolff
-cluster and, since the detailed balance with respect to the spin--spin energy
-is performed by the cluster creation, only the change in energy due to the
-external field is considered. For $H=0$ this reverts to the Wolff algorithm
-exactly. A strictly faster version of this algorithm is accomplished by
-evaluating the probability that the cluster-building is reverted at every step
-in its formation, so that large unlikely clusters abandoned before very much
-time is wasted in constructing them. This algorithm is given as follows.
+The Wolff algorithm 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^{-\beta
+H}$ every time a spin is added to it.
 
 \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$ be the index of a random site
-      \STATE $q.\mathrm{push}(i_0)$
+      \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.\mathrm{pop}$
+        \STATE $i\gets q.\mathtt{pop}$
         \IF {$s_i=\sigma$}
-          \FOR {$j=1$ \TO $n$}
-            \IF {$1-e^{-2\beta J_{ij}s_is_j}\geq$ random number
-          between 0 and 1}
-          \STATE $q.\mathrm{push}(j)$
+          \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$
@@ -239,35 +223,36 @@ time is wasted in constructing them. This algorithm is given as follows.
 \end{figure}
 
 
+
 \begin{figure}
   \begin{algorithm}[H]
     \begin{algorithmic}
       \REQUIRE $s\in\G^n$
       \STATE \textbf{let} $q_1$ and $q_2$ be empty stacks
-      \STATE \textbf{let} $i_0$ be the index of a random site
-      \STATE $q_1.\mathrm{push}(i_0)$
+      \STATE \textbf{let} $i_0\in V$
+      \STATE $q_1.\mathtt{push}(i_0)$
       \STATE $\sigma\gets s_{i_0}$
-      \STATE stopped $\gets$ \FALSE
+      \STATE \texttt{stopped} $\gets$ \textbf{false}
       \WHILE {$q_1$ is not empty}
-        \STATE $i\gets q_1.\mathrm{pop}$
+        \STATE $i\gets q_1.\mathtt{pop}$
         \IF {$s_i=\sigma$}
-          \FORALL {$j>0$ such that $J_{ij}\neq0$}
-            \IF {$s_is_j=1$ \AND $1-e^{-2\beta J_{ij}}\geq$ random number
-          between 0 and 1}
-          \STATE $q_1.\mathrm{push}(j)$
+          \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_1.\mathtt{push}(j)$
             \ENDIF
           \ENDFOR
           \STATE $s_i\gets-s_i$
-          \STATE $q_2.\mathrm{push}(i)$
-          \IF {$1-e^{-2\beta\sigma H}>$ random number between 0 and 1}
-            \STATE stopped $\gets$ \TRUE
+          \STATE $q_2.\mathtt{push}(i)$
+          \IF {$1-\exp(-2\beta\sigma H)>\mathop{\mathtt{UniformRandom}}(0,1)$}
+          \STATE \texttt{stopped} $\gets$ \textbf{true}
             \STATE \textbf{break}
           \ENDIF
         \ENDIF
       \ENDWHILE
       \IF {stopped}
         \WHILE {$q_2$ is not empty}
-          \STATE $i\gets q_2.\mathrm{pop}$
+          \STATE $i\gets q_2.\mathtt{pop}$
           \STATE $s_i\gets-s_i$
         \ENDWHILE
       \ENDIF
@@ -279,102 +264,161 @@ time is wasted in constructing them. This algorithm is given as follows.
 
 \section{Cluster-flip in a Field}
 
-Consider a system whose state $s=\{s_0,s_1,\ldots,s_n\}$ $s\in\Z_2^{n+1}$ is specified by the
-value of $n+1$ boolean variables with the Hamiltonian
+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. 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
 \[
-  \new\H(s)=-\sum_{i=0}^n\sum_{j=0}^n\new J_{ij}s_is_j,
-  \label{eq:new.ham}
+  \tilde J_{ij}=\begin{cases}
+    J_{ij} & i,j>0\\
+    H & \text{otherwise},
+  \end{cases}
 \]
-where
+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$. Now define a
+Hamiltonian function 
+$\tilde\H:S^{n+1}\to\R$ on this new, larger configuration space defined by
 \[
-  \new J_{ij}=\begin{cases}
-    1 & \text{$i,j>0$, $i$, $j$ are neighboring sites}\\
-    H & \text{$i=0$ or $j=0$}\\
-    0 & \text{otherwise}
-  \end{cases}
+  \tilde\H(s)=-\sum_{\{i,j\}\in\tilde E}\tilde J_{ij}s_is_j.
 \]
-gives the coupling. Since \eqref{eq:new.ham} is a discrete-spin Hamiltonian
-invariant under the transformation $s\to-s$, it directly admits analysis using
-the Wolff algorithm specified above. By rewriting the sums, however, we also
-see that
+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 write
 \[
-  \new\H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-H\new M(s)
-  \label{eq:mod_hamiltonian}
+  \tilde\H(s)=-\sum_{\{i,j\}\in E}J_{ij}s_is_j-H\tilde M(s)
 \]
-where
+where the new magnetization $\tilde M:S^{n+1}\to\R$ is defined for $(s_0,s)\in
+S^{n+1}$ by
 \[
-  \new M(s)=s_0\sum_{i=1}^ns_i
+  \tilde M(s)=s_0\sum_{i\in V}s_i=M(s_0\times(s_1,\ldots,s_n)).
 \]
-This is precisely the same Hamiltonian as \eqref{eq:ham} but with a redefined
-magnetization whose sign depends on the state of the zeroth spin. Our
-algorithm for sampling the Ising model in a field, then, is simply the Wolff
-algorithm described above with $\tilde J_{ij}$ replacing $J_{ij}$. These are
-formally equivalent, since every state $s\in S$ corresponds to two states
-$\tilde s\in\tilde S$, $\tilde s=\pm(1,s)$, such that $\H(s)=\tilde\H(\tilde
-s)$. We therefore have
-\begin{align}
-  \avg A_{\tilde S}
-  &=\frac{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde
-  s)}A(s_0s)}{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde s)}}\\
-  &=\frac{\sum_se^{-\beta\tilde\H(s_0,s)}A(s)+\sum_se^{-\beta\tilde\H(-s_0,-s)}A(s)}{\sum_se^{-\beta\tilde\H(1,s)}+\sum_se^{-\beta\tilde\H(-1,-s)}}\\
-  &=\avg A_S
-\end{align}
-$\tilde A(s_0,s)=A(s_0s)$. This is true for the mangetization $\tilde M$ as
-defined above, for instance.
+In fact, any observable $A$ of the original system can be written as an
+observable $\tilde A$ of the new system by defining
+\[
+  \tilde A(s)=A(s_0\times(s_1,\ldots,s_n))
+\]
+and the expected value of the observable $A$ in the old system and that of its
+counterpart $\tilde A$ in the new system is unchanged. This can be seen by
+using the facts that
+$\tilde\H(s)=\tilde\H(-s)$, $\sum_{s\in S^n}f(s)=\sum_{s\in S^n}F(-s)$ for any
+$f:S^n\to\R$, and $\tilde\H(1,s)=\H(s)$ for $s\in S^n$, from which follows
+that
+\[
+  \begin{aligned}
+  \tilde Z\avg{\tilde A}
+  &=\sum_{s\in S^{n+1}}\tilde A(s)e^{-\beta\tilde\H(s)}
+  =\sum_{s_0\in S}\sum_{s\in S^n}\tilde A(s_0,s)e^{-\beta\tilde\H(s_0,s)}\\
+  &=\bigg(\sum_{s\in S^n}A(s)e^{-\beta\tilde\H(1,s)}
+    +\sum_{s\in S^n}A(-s)e^{-\beta\tilde\H(-1,s)}\bigg)\\
+  &=\bigg(\sum_{s\in S^n}A(s)e^{-\beta\tilde\H(1,s)}
+    +\sum_{s\in S^n}A(-s)e^{-\beta\tilde\H(1,-s)}\bigg)\\
+  &=\bigg(\sum_{s\in S^n}A(s)e^{-\beta\tilde\H(1,s)}
+    +\sum_{s\in S^n}A(s)e^{-\beta\tilde\H(1,s)}\bigg)\\
+  &=\bigg(\sum_{s\in S^n}A(s)e^{-\beta\H(s)}
+    +\sum_{s\in S^n}A(s)e^{-\beta\H(s)}\bigg)\\
+    &=2Z\avg A.
+\end{aligned}
+\]
+An identical calculation shows $\tilde Z=2Z$, therefore immediately proving
+$\avg{\tilde A}=\avg A$. Notice this correspondence also holds for the
+Hamiltonian.
+
+Our new spin system with an additional field is, when $H$ is greater than
+zero, 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}.
 
 \begin{figure}
   \begin{algorithm}[H]
     \begin{algorithmic}
-      \REQUIRE $s\in\G^n$
+      \REQUIRE $s\in S^n$
       \STATE \textbf{let} $q$ be an empty stack
-      \STATE \textbf{let} $i_0$ be the index of a random site
-      \STATE $q.\mathrm{push}(i_0)$
+      \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.\mathrm{pop}$
+        \STATE $i\gets q.\mathtt{pop}$
         \IF {$s_i=\sigma$}
-          \FOR {$j=0$ \TO $n$}
-            \IF {$1-\exp(-2\beta \tilde J_{ij}s_is_j)\geq$ random number
-          between 0 and 1}
-          \STATE $q.\mathrm{push}(j)$
+          \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 (Non-zero Field)}
+    \caption{Wolff (Nonzero Field)}
     \label{alg:wolff-field}
   \end{algorithm}
 \end{figure}
 
 
 At sufficiently small $H$ both our algorithm and the hybrid Wolff--Metropolis reduce
-to the ordinary Wolff algorithm. At very large $H$, the hybrid
+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 the field.
-
-We measured the correlation time $\tau$ of these three algorithms at various fields
-and temperatures. This was done using a batch mean. Time was measured as ``the
-number of spins that have either been attempted to flip or attempted to add to
-a cluster.''
-
-$E_s=A_s+XB_s$
-\[
-  \pd{\avg B}X=\beta\big(\avg{B^2}-\avg B^2\big)
-\]
-\[
-  \chi=\pd{\avg M}H=\beta\big(\avg{M^2}-\avg M^2\big)
-\]
-
-Since \eqref{eq:mod_hamiltonian} has the form of a simple (ferromagnetic)
-spin-spin Hamiltonian, its ensemble average can be taken using the Wolff
-algorithm while never rejecting a step.
+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 a batch mean estimator. 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}
+\end{figure}
 
 \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 
+
 For any finite system, the average magnetization at zero field is identically
 zero at all temperatures. However, people often want to use finite simulations
 to estimate the average magnetization of the Ising model in the thermodynamic