many small revisions

2 files changed, 198 insertions, 134 deletions

diff --git a/monte-carlo.bib b/monte-carlo.bib
index 0ca69fb..6eb96d1 100644
--- a/monte-carlo.bib
+++ b/monte-carlo.bib
@@ -411,6 +411,22 @@
 	file = {APS Snapshot:/home/pants/.zotero/data/storage/Q5RL6Q8U/PhysRevE.58.html:text/html;Redner et al. - 1998 - Graphical representations and cluster algorithms f.pdf:/home/pants/.zotero/data/storage/WYEU9G6Y/Redner et al. - 1998 - Graphical representations and cluster algorithms f.pdf:application/pdf}
 }
 
+@article{wolff_lattice_1988,
+	title = {Lattice field theory as a percolation process},
+	volume = {60},
+	url = {https://link.aps.org/doi/10.1103/PhysRevLett.60.1461},
+	doi = {10.1103/PhysRevLett.60.1461},
+	abstract = {For a given lattice spin or gauge theory an associated correlated bond or plaquette percolation process is constructed. It is conjectured to reproduce the universal scaling behavior of the original model. Different field theories lead to different cluster weights generalizing a result by Fortuin and Kasteleyn for Potts models. The new representation lends itself to the design of Monte Carlo algorithms with reduced critical slowing down.},
+	number = {15},
+	urldate = {2018-04-17},
+	journal = {Physical Review Letters},
+	author = {Wolff, Ulli},
+	month = apr,
+	year = {1988},
+	pages = {1461--1463},
+	file = {APS Snapshot:/home/pants/.zotero/data/storage/EFPZ9VED/PhysRevLett.60.html:text/html;Wolff - 1988 - Lattice field theory as a percolation process.pdf:/home/pants/.zotero/data/storage/EBMJTQ4T/Wolff - 1988 - Lattice field theory as a percolation process.pdf:application/pdf}
+}
+
 @article{guida_critical_1998,
 	title = {Critical exponents of the {N} -vector model},
 	volume = {31},
@@ -428,6 +444,24 @@
 	file = {IOP Full Text PDF:/home/pants/.zotero/data/storage/K468APXL/Guida and Zinn-Justin - 1998 - Critical exponents of the N -vector model.pdf:application/pdf}
 }
 
+@article{fortuin_random-cluster_1972,
+	title = {On the random-cluster model: {I}. {Introduction} and relation to other models},
+	volume = {57},
+	issn = {0031-8914},
+	shorttitle = {On the random-cluster model},
+	url = {http://www.sciencedirect.com/science/article/pii/0031891472900456},
+	doi = {10.1016/0031-8914(72)90045-6},
+	abstract = {The random-cluster model is defined as a model for phase transitions and other phenomena in lattice systems, or more generally in systems with a graph structure. The model is characterized by a (probability) measure on a graph and a real parameter κ. By specifying the value of κ to 1, 2, 3, 4, … is shown that the model covers the percolation model, the Ising model, the Ashkin-Teller-Potts model with 3, 4, … states per atom, respectively, and thereby, contains information on graph-colouring problems; in the limit κ ↓ 0 it describes linear resistance networks. It is shown that the function which for the random-cluster model plays the role of a partition function, is a generalization of the dichromatic polynomial earlier introduced by Tutte, and related polynomials.},
+	number = {4},
+	urldate = {2018-04-20},
+	journal = {Physica},
+	author = {Fortuin, C. M. and Kasteleyn, P. W.},
+	month = feb,
+	year = {1972},
+	pages = {536--564},
+	file = {ScienceDirect Full Text PDF:/home/pants/.zotero/data/storage/463W5CMX/Fortuin and Kasteleyn - 1972 - On the random-cluster model I. Introduction and r.pdf:application/pdf;ScienceDirect Snapshot:/home/pants/.zotero/data/storage/EEHFQXUY/0031891472900456.html:text/html}
+}
+
 @article{stauffer_scaling_1979,
 	title = {Scaling theory of percolation clusters},
 	volume = {54},
diff --git a/monte-carlo.tex b/monte-carlo.tex
index fe6227d..9b61ba1 100644
--- a/monte-carlo.tex
+++ b/monte-carlo.tex
@@ -53,13 +53,15 @@
 
 \begin{abstract}
   We introduce a natural way to extend celebrated spin-cluster Monte Carlo
-  algorithms for fast thermal lattice simulations at criticality, like Wolff, to
-  systems in arbitrary fields. The method relies on the generalization of the
-  `ghost spin' representation to one with a `ghost transformation' that
-  restores invariance to spin symmetries at the cost of an extra degree of
+  algorithms for fast thermal lattice simulations at criticality, like Wolff,
+  to systems in arbitrary fields. By generalizing the `ghost spin'
+  representation to one with a `ghost transformation,\!' global invariance to
+  spin symmetry transformations is restored at the cost of an extra degree of
   freedom. The ordinary cluster-building process can then be run on the new
-  representation. For several canonical systems, we show that this extension
-  preserves the scaling of accelerated dynamics in the absence of a field.
+  representation. We show that this extension preserves the scaling of
+  accelerated dynamics in the absence of a field for several canonical systems
+  and demonstrate the method's use in modelling the presence of novel
+  nonlinear fields.
 \end{abstract}
 
 \maketitle
@@ -68,9 +70,9 @@ Lattice models are important in the study of statistical physics and phase
 transitions. Rarely exactly solvable, they are typically studied by
 approximate and numerical methods. Monte Carlo techniques are a common way of
 doing this, approximating thermodynamic quantities by sampling the
-distribution of systems states. These Monte Carlo algorithms are better the
-faster they arrive at a statistically independent sample. This typically
-becomes a problem near critical points, where critical slowing down
+distribution of system states. These Monte Carlo algorithms are better the
+faster they arrive at a statistically independent sample. This becomes a
+problem near critical points, where critical slowing down
 \cite{wolff_critical_1990} results in power-law divergences of dynamic
 timescales.
 
@@ -82,8 +84,8 @@ dynamic exponents for many spin systems \cite{wolff_comparison_1989,
 du_dynamic_2006, liu_dynamic_2014, wang_cluster_1990}, including the Ising,
 $\mathrm O(n)$ \cite{wolff_collective_1989}, and Potts
 \cite{swendsen_nonuniversal_1987, baillie_comparison_1991} models. These
-algorithms rely on the natural symmetry of the systems in question under
-symmetry operations on the spins.
+algorithms rely on the natural invariance of the systems in question under
+symmetry transformations on their spins.
 
 Some success has been made in extending these algorithms to systems in certain
 external fields by adding a `ghost site' \cite{coniglio_exact_1989} that
@@ -97,7 +99,7 @@ ala-nissila_numerical_1994}, and other categories of fields have been applied
 using replica methods \cite{redner_graphical_1998, chayes_graphical_1998,
 machta_replica-exchange_2000}.  Monte Carlo techniques that involve cluster
 updates at fixed magnetization have been used to examine quantities at fixed
-field by integrating the associated thermodynamic functions
+field by later integrating measured thermodynamic functions
 \cite{martin-mayor_cluster_2009, martin-mayor_tethered_2011}.
 
 We show that the scaling of correlation time near the critical point of
@@ -107,28 +109,32 @@ field to various non-symmetric perturbations. We also show, by a redefinition
 of the spin--spin coupling in a generic class of spin systems,
 \emph{arbitrary} external fields can be treated using cluster methods. Rather
 than the introduction of a `ghost spin,\!' our representation relies on
-introducing a `ghost transformation.\!'
+introducing a `ghost transformation,\!' an extra degree of freedom residing on
+a `ghost' site coupled to all other sites that takes its values from the
+collection of spin symmetry transformations of the base model rather than
+resemble the base spins themselves.
 
 \section{Clusters Without a Field}
 
 We will pose the problem in a general way, but several specific examples can
-be found in Table~\ref{table:models} for concreteness. Let $G=(V,E)$ be a graph, where the set of vertices $V=\{1,\ldots,N\}$
-enumerates the sites of a lattice and the set of edges $E$ contains pairs of
-neighboring sites. Let $R$ be a group acting on a set $X$, with the action of
-group elements $r\in R$ on elements $s\in X$ denoted $r\cdot s$. $X$ is the
-set of states accessible by a spin, and $R$ is the \emph{symmetry group} of
-$X$. The set $X$ must admit a measure $\mu$ that is invariant under the action
-of $R$, e.g., for any $A\subseteq X$ and $r\in R$, $\mu(r\cdot A)=\mu(A)$.
-This trait is shared by the counting measure on any discrete set, or by any
-group acting by isometries on a Riemannian manifold, such as $\mathrm O(n)$ on
-$S^{n-1}$ in the $\mathrm O(n)$ model \cite{caracciolo_wolff-type_1993}.
-Finally, a subset $R_2$ of elements in $R$ of order two must act transitively on
-$X$. This property, while apparently obscure, is shared by any symmetric space
+be found in Table~\ref{table:models} for concreteness. Let $G=(V,E)$ be a
+graph, where the set of vertices $V=\{1,\ldots,N\}$ enumerates the sites of a
+lattice and the set of edges $E$ contains pairs of neighboring sites. Let $R$
+be a group acting on a set $X$, with the action of group elements $r\in R$ on
+elements $s\in X$ denoted $r\cdot s$. $X$ is the set of states accessible by
+each 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)$ models \cite{caracciolo_wolff-type_1993}.  Finally, a subset $R_2$ 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
+group. In fact, all the examples listed here have spin spaces with natural
 metrics whose symmetry group is their set of isometries.  We put one spin at
 each site of the lattice described by $G$, so that the state of the entire
-spin system is described by elements $\set s\in X\times\cdots\times X=X^N$. 
+system is described by elements $\set s\in X\times\cdots\times X=X^N$. 
 
 The Hamiltonian of this system is a function $\H:X^N\to\R$ defined by
 \[
@@ -143,7 +149,7 @@ interactions---or allow $B$ to be a function of site---for applying arbitrary
 boundary conditions or modelling random fields. The formal results of this
 paper (that the algorithm obeys detailed balance and ergodicity) hold equally
 well for these cases, but we will drop the additional index notation for
-clarity. Statements about efficiency may not.
+clarity. Statements about efficiency may not hold.
 
 \begin{table*}[htpb]
   \begin{tabular}{l||ccccc}
@@ -200,12 +206,12 @@ in the following way.
 \begin{enumerate}
   \item Pick a random site $m_0$ and add it to the stack.
 
-  \item Select a rotation $r\in R_2$ distributed by $f(r\mid m_0,\set s)$.
+  \item Select a transformation $r\in R_2$ distributed by $f(r\mid m_0,\set s)$.
     Often $f$ is taken as uniform on $R_2$, but it is sufficient for preserving
     detailed balance that $f$ be any function of the seed site $m_0$ and
     $Z(s,r\cdot s)$ for all $s\in\set s$. The flexibility offered by the
-    choice of distribution will be useful in situations where the state space
-    is infinite.
+    choice of distribution will be useful in situations where the set of spin
+    states is infinite.
   \item While the stack isn't empty,
     \begin{enumerate}
       \item pop site $m$ from the stack.
@@ -223,16 +229,16 @@ in the following way.
     \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
+transformed by the action of $r$. In order for this algorithm to be useful, it
 must satisfy ergodicity and detailed balance. Ergodicity is satisfied since we
-have ensured that the $R_2$ acts transitively on $X$, e.g., for any $s,t\in X$
+have ensured that $R_2$ acts transitively on $X$, e.g., for any $s,t\in X$
 there exists $r\in R_2$ such that $r\cdot s=t$. Since there is a nonzero
-probability that only one spin is rotated and that spin can be rotated into
-any state, ergodicity follows. The probability $P(\set s\to\set{s'})$ that the
-configuration $\set s$ is brought to $\set 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(\set{s'}\to\set s)$ by
+probability that only one spin is transformed and that spin can be transformed
+into any state, ergodicity follows. The probability $P(\set s\to\set{s'})$
+that the configuration $\set s$ is brought to $\set s'$ by the flipping of a
+cluster formed by accepting transformations of spins via bonds $C\subseteq E$
+and rejecting transformations via bonds $\partial C\subset E$ is related to
+the probability of the reverse process $P(\set{s'}\to\set s)$ by
 \begin{widetext}
 \[
   \begin{aligned}
@@ -250,7 +256,7 @@ s',s')=Z(r\cdot s,s)$.
 
 The Wolff algorithm is well known to be efficient in sampling many spin models
 near and away from criticality, including the Ising, Potts, and $\mathrm O(n)$
-models. In general, its efficiently will depend on the system at hand, e.g.,
+models. In general, its efficiency will depend on the system at hand, e.g.,
 the structure of the configurations $X$ and group $R$. A detailed discussion
 of this dependence for a class of configuration spaces with continuous
 symmetry groups can be found in \cite{caracciolo_generalized_1991,
@@ -299,8 +305,8 @@ $r,s_0\in R$ and $s\in X$,
   =\tilde\J(s_0,s)
 \end{aligned}
 \]
-The invariance of $\tilde\J$ to rotations given other arguments follows from
-the invariance properties of $\J$.
+The invariance of $\tilde\J$ to global transformations given other arguments
+follows from the invariance properties of $\J$.
 
 We have produced a system incorporating the field function $B$ whose
 Hamiltonian is invariant under global rotations, but how does it relate to our
@@ -358,6 +364,24 @@ value of $\tilde A$ on the new system. In contrast with the simpler ghost spin
 representation, this form of the Hamiltonian might be considered the `ghost
 transformation' representation.
 
+One of the celebrated features of the cluster representation of the Ising and
+associated models are the improved estimators of various quantities in the
+base model, found by measuring conjugate properties of the clusters themselves
+\cite{wolff_lattice_1988}. What of these quantities survive this translation?
+As is noted in the formative construction of the cluster representation for
+the Ising and Potts models, all estimators involving correlators between spins
+are preserved, including correlators with the ghost site
+\cite{fortuin_random-cluster_1972}. Where a previous improved estimator
+exists, we expect this representation to extend it to finite field, all other
+features of the algorithm held constant. For instance, the average cluster
+size in the Wolff algorithm is often said to be an estimator for the magnetic
+susceptibility in the Ising and Potts models (average sum over $x$-components
+in the cluster for the $\mathrm O(n)$), but really what it estimates is the
+averaged squared magnetization, which corresponds to the susceptibility when
+the average magnetization is zero. At finite field the latter thing is no
+longer true, but the correspondence between cluster size (sum of
+$x$-components) and the squared magnetization continues to hold.
+
 \section{Examples}
 \label{sec:examples}
 
@@ -385,59 +409,61 @@ authors in an existing interactive Ising simulator at
 \label{sec:examples:on}
 
 In the $\mathrm O(n)$ model spins are described by vectors on the
-$(n-1)$-sphere $S^{n-1}$. Its symmetry group is $O(n)$, $n\times n$ orthogonal
-matrices, which act on the spins by matrix multiplication. The elements of
-$O(n)$ of order two are reflections about hyperplanes through the origin and
-$\pi$ rotations about any axis through the origin. Since the former generate
-the entire group, reflections alone suffice to provide ergodicity. Sampling
-those reflections uniformly works well at criticality. The `ghost spin'
-version of the algorithm has been used to apply a simple vector field to the
-$\mathrm O(3)$ model \cite{dimitrovic_finite-size_1991}. Other fields of
-interest include $(n+1)$-dimensional spherical harmonics
-\cite{jose_renormalization_1977} and cubic fields \cite{bruce_coupled_1975,
-blankschtein_fluctuation-induced_1982}, which can be applied with the new
-method. The method is quickly generalized to spins whose symmetry groups are
-other compact Lie groups \cite{caracciolo_generalized_1991,
-caracciolo_wolff-type_1993}.
-
-At low temperature or high external vector field field selecting reflections
+$(n-1)$-sphere $S^{n-1}$. Its symmetry group is $\mathrm O(n)$, $n\times n$
+orthogonal matrices, which act on the spins by matrix multiplication. The
+elements of $\mathrm O(n)$ of order two are reflections about hyperplanes
+through the origin and $\pi$ rotations about any axis through the origin.
+Since the former generate the entire group, reflections alone suffice to
+provide ergodicity. Sampling those reflections uniformly works well at
+criticality. The `ghost spin' version of the algorithm has been used to apply
+a simple vector field to the $\mathrm O(3)$ model
+\cite{dimitrovic_finite-size_1991}. Other fields of interest include
+$(n+1)$-dimensional spherical harmonics \cite{jose_renormalization_1977} and
+cubic fields \cite{bruce_coupled_1975, blankschtein_fluctuation-induced_1982},
+which can be applied with the new method. The method is quickly generalized to
+spins whose symmetry groups are other compact Lie groups
+\cite{caracciolo_generalized_1991, caracciolo_wolff-type_1993}.
+
+At low temperature or high external vector field selecting reflections
 uniformly becomes inefficient because the excitations of the model are spin
 waves, in which the magnetization only differs by a small amount between
 neighboring spins. Under these conditions, most choices of reflection plane
 will cause a change in energy so great that the whole system is always
 flipped, resulting in many correlated samples. To ameliorate this, one can
 draw reflections from a distribution that depends on how the seed spin is
-transformed. We implement this in the following way. Say that the seed of the
-cluster is $s$. Generate a vector $t$ taken uniformly from the space of unit
-vectors orthogonal to $s$. Let the plane of reflection be that whose normal is
-$n=s+\zeta t$, where $\zeta$ is drawn from a normal distribution of mean zero
-and variance $\sigma$. It follows that the tangent of the angle between $s$
-and the plane of reflection is also distributed normally with zero mean and
-variance $\sigma$. Since the distribution of reflection planes only depends on
-the angle between $s$ and the plane, and since that angle is invariant under
-the reflection, this choice preserves detailed balance.
+transformed. We implement this in the following way. Say that the state of the
+seed of the cluster is $s$. Generate a vector $t$ taken uniformly from the
+space of unit vectors orthogonal to $s$. Let the plane of reflection be that
+whose normal is $n=s+\zeta t$, where $\zeta$ is drawn from a normal
+distribution of mean zero and variance $\sigma$. It follows that the tangent
+of the angle between $s$ and the plane of reflection is also distributed
+normally with zero mean and variance $\sigma$. Since the distribution of
+reflection planes only depends on the angle between $s$ and the plane, and
+since that angle is invariant under the reflection, this choice preserves
+detailed balance.
 
 The choice of $\sigma$ can be inspired by mean field theory. At high field or
 low temperature, spins are likely to both align with the field and each other
-and the model is asymptotically equal to a simple Gaussian one, with in the
-limit of large $L$ the expected square angle between neighbors being
+and the model is asymptotically equal to a simple Gaussian one, in which in the
+limit of large $L$ the expected square angle between neighbors is
 \[
   \avg{\theta^2}\simeq\frac{(n-1)T}{D+H/2}.
 \]
 We take $\sigma=\sqrt{\avg{\theta^2}}/2$. Fig.~\ref{fig:generator_times} shows
-the effect of making such a choice on autocorrelation times for a critical
-\threedee \textsc{xy} ($\mathrm O(2)$) model. At small fields both methods
-perform the same as zero field Wolff.  Intermediate field values see
+the effect of making such a choice on autocorrelation times for the energy for
+a critical \threedee \textsc{xy} ($\mathrm O(2)$) model. At small fields both
+methods perform the same as zero field Wolff.  Intermediate field values see
 efficiency gains for both methods. At large field the uniform sampling method
-sees correlation times grow rapidly, while for the sampling method described
-here the correlation time crosses over to a constant. A similar behavior holds
-for the critical $\mathrm O(3)$ model, though in that case the constant value
-the correlation time approaches at large field is larger than its minimum
-value (see Fig.~\ref{fig:correlation_time-collapse}). This behavior isn't
-particularly worrisome, since the very large field regime corresponds to
-correlation lengths smaller than the lattice spacing and is well-described by
-other algorithms.  More detailed discussion on correlation times and these
-numeric experiments can be found in section \ref{sec:performance}.
+sees correlation times grow rapidly without bound, while for the sampling
+method described here the correlation time crosses over to a constant. A
+similar behavior holds for the critical $\mathrm O(3)$ model, though in that
+case the constant value the correlation time approaches at large field is
+larger than its minimum value (see Fig.~\ref{fig:correlation_time-collapse}).
+This behavior isn't particularly worrisome, since the very large field regime
+corresponds to correlation lengths smaller than the lattice spacing and is
+well-described by other algorithms.  More detailed discussion on correlation
+times and these numeric experiments can be found in section
+\ref{sec:performance}.
 
 \begin{figure}
   \include{fig_generator-times}
@@ -485,13 +511,13 @@ symmetry group is the infinite dihedral group $D_\infty=\{r_i,s_i\mid
 i\in\Z\}$, whose action on the spin $j\in\Z$ is given by $r_i\cdot j=i+j$ and
 $s_i\cdot j=-i-j$. The elements of order two are reflections $s_i$, whose
 action on $\Z$ is transitive. The coupling can be any function of the absolute
-difference $|i-j|$.  Because random choice of reflection will almost always
+difference $|i-j|$.  Because uniform choice of reflection will almost always
 result in energy changes so large that the whole system is flipped, it is
 better to select random reflections about integers or half-integers close to
-the average state of the system.  A variant of the algorithm has been applied
-without a field whose success relies both on this and another technique
+the state of the system.  A variant of the algorithm has been applied without
+a field whose success relies both on this and another technique
 \cite{evertz_stochastic_1991}. They note that detailed balance is still
-satisfied if the bond probabilities \eqref{eq:bond_probability} is modified by
+satisfied if the bond probabilities \eqref{eq:bond_probability} are modified by
 adding a constant $0<x\leq1$ with
 \[
   p_r(s_m,s_j\mid x)=\min\{0,1-xe^{\beta(\J(r\cdot s_m,s_j)-\J(s_m,s_j))}\}.
@@ -530,31 +556,6 @@ one should expect its natural extension in the presence of a field to scale
 roughly like $h^{-z\nu/\beta\delta}$ and collapse appropriately as a function
 of $hL^{\beta\delta/\nu}$.
 
-We measured the autocorrelation time $\tau$ of the energy $\H$ for a variety of
-models at critical temperature with many system sizes and canonical fields
-(see Table~\ref{table:models} with $h=\beta H$) using standard methods for
-obtaining the value and uncertainty from timeseries
-\cite{ossola_dynamic_2004}. Since the computational effort expended in each
-step of the algorithm depends linearly on the size of the associated cluster,
-these values are then scaled by the average cluster size per site
-$\avg{s_{\text{\sc 1c}}}/L^D$ to produce something proportional to machine
-time. The resulting scaling behavior, plotted in
-Fig.~\ref{fig:correlation_time-collapse}, is indeed consistent with an
-extension to finite field of the behavior at zero field, with an eventual
-finite-size crossover to constant autocorrelation time. This crossover isn't
-always kind to the efficiency, e.g., in the $\mathrm O(3)$ model, but in the
-large-field regime where the crossover happens the correlation length is on
-the scale of the lattice spacing and better algorithms exist, like
-Bortz--Kalos--Lebowitz for the Ising model \cite{bortz_new_1975}. Also plotted
-are lines proportional to $h^{-z\nu/\beta\delta}$, which match the behavior of
-the correlation times in the intermediate scaling region. Values of the
-critical exponents for the models were taken from the literature
-\cite{wu_potts_1982, el-showk_solving_2014, guida_critical_1998} with the
-exception of $z$ for the energy in the Wolff algorithm, which was determined
-for each model by making a power law fit to the constant low field behavior.
-These exponents are imprecise and are provided with only qualitative
-uncertainty.
-
 \begin{figure*}
   \include{fig_correlation-times}
   \caption{
@@ -566,19 +567,48 @@ uncertainty.
     $h^{-z\nu/\beta\delta}$ for each model. The dynamic exponents $z$ are
     roughly measured as \twodee Ising: 0.23(2), \threedee Ising: 0.28(2),
     \twodee 3-State Potts: 0.55(1), \twodee 4-State Potts: 0.94(5),
-    \threedee O(2): 0.17(2), \threedee O(3): 0.13(2).
+    \threedee O(2): 0.17(2), \threedee O(3): 0.13(2). $\mathrm O(n)$ models
+    use the distribution of transformations described in Section
+    \ref{sec:examples:on}.
   }
   \label{fig:correlation_time-collapse}
 \end{figure*}
 
+We measured the autocorrelation time $\tau$ of the energy $\H$ for a variety
+of models at critical temperature with many system sizes and canonical fields
+(see Table~\ref{table:models} with $h=\beta H$) using standard methods for
+obtaining the value and uncertainty from timeseries
+\cite{ossola_dynamic_2004}. Since the computational effort expended in each
+step of the algorithm depends linearly on the size of the associated cluster,
+these values are then scaled by the average cluster size per site
+$\avg{s_{\text{\sc 1c}}}/L^D$ to produce something proportional to machine
+time per site. The resulting scaling behavior, plotted in
+Fig.~\ref{fig:correlation_time-collapse}, is indeed consistent with an
+extension to finite field of the behavior at zero field, with an eventual
+finite-size crossover to constant autocorrelation time at large field. This
+crossover isn't always kind to the efficiency, e.g., in the $\mathrm O(3)$
+model, but in the large-field regime where the crossover happens the
+correlation length is on the scale of the lattice spacing and better
+algorithms exist, like Bortz--Kalos--Lebowitz for the Ising model
+\cite{bortz_new_1975}. Also plotted are lines proportional to
+$h^{-z\nu/\beta\delta}$, which match the behavior of the correlation times in
+the intermediate scaling region as expected. Values of the critical exponents
+for the models were taken from the literature \cite{wu_potts_1982,
+el-showk_solving_2014, guida_critical_1998} with the exception of $z$ for the
+energy in the Wolff algorithm, which was determined for each model by making a
+power law fit to the constant low field behavior.  These exponents are
+imprecise and are provided in the figure with only qualitative uncertainty.
+
 Since the formation and flipping of clusters is the hallmark of Wolff
 dynamics, another way to ensure that the dynamics with field scale like those
 without is to analyze the distribution of cluster sizes. The success of the
 algorithm at zero field is related to the fact that the clusters formed
 undergo a percolation transition at models' critical point.  According to the
 scaling theory of percolation \cite{stauffer_scaling_1979}, the distribution
-of cluster sizes in a full Swendsen--Wang decomposition of the system scales
-consistently near the critical point if it has the form
+of cluster sizes in a full Swendsen--Wang decomposition---where the whole
+system is decomposed into clusters with every bond activated with probability
+\eqref{eq:bond_probability}---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}).
 \]
@@ -593,6 +623,18 @@ proportional to their size, or
   \end{aligned}
 \]
 
+\begin{figure*}
+  \input{fig_clusters_ising2d}
+  \caption{Collapses of rescaled average Wolff cluster size $\avg s_{\text{\sc
+    1c}}L^{-\gamma/\nu}$ as a function of field scaling variable
+    $hL^{\beta\delta/\nu}$ for a variety of models. Critical exponents
+    $\gamma$, $\nu$, $\beta$, and $\delta$ are model-dependant. Colored lines
+    and points depict values as measured by the extended algorithm. Solid
+    black lines show a plot of $g(0,x)\propto x^{2/\delta}$ for each model.
+  }
+  \label{fig:cluster_scaling}
+\end{figure*}
+
 For the Ising model, an additional scaling relation can be written. Since the
 average cluster size is the average squared magnetization, it can be related
 to the scaling functions of the magnetization and susceptibility per site by
@@ -619,35 +661,23 @@ cluster size collapses for each model according to the scaling hypothesis, and
 the large-field behavior likewise scales as we expect from the na\"ive Ising
 conjecture.
 
-\begin{figure*}
-  \input{fig_clusters_ising2d}
-  \caption{Collapses of rescaled average Wolff cluster size $\avg s_{\text{\sc
-    1c}}L^{-\gamma/\nu}$ as a function of field scaling variable
-    $hL^{\beta\delta/\nu}$ for a variety of models. Critical exponents
-    $\gamma$, $\nu$, $\beta$, and $\delta$ are model-dependant. Colored lines
-    and points depict values as measured by the extended algorithm. Solid
-    black lines show a plot of $g(0,x)\propto x^{2/\delta}$ for each model.
-  }
-  \label{fig:cluster_scaling}
-\end{figure*}
-
 \section{Applying Nonlinear Fields to the xy Model}
 
 Thus far our numeric work has quantified the performance of existing
 techniques. Briefly, we demonstrate our general framework in a new way:
 harmonic perturbations to the low-temperature {\sc xy}, or \twodee O(2),
 model. We consider fields of the form $B_n(s)=h_n\cos(n\theta(s))$, where
-$\theta$ is the angle made between $s$ and the $x$-axis, say. Corrections of
-these types are expected to appear in realistic models of systems na\"ively
-expected to exhibit Kosterlitz--Thouless critical behavior due to the presence
-of the lattice or substrate. Whether these fields are relevant or irrelevant
-in the renormalization group sense determines whether those systems spoil or
-admit that critical behaviour. Among many fascinating
+$\theta$ is the angle made between $s$ and the $x$-axis. Corrections of these
+types are expected to appear in realistic models of systems na\"ively expected
+to exhibit Kosterlitz--Thouless critical behavior due to the presence of the
+lattice or substrate. Whether these fields are relevant or irrelevant in the
+renormalization group sense determines whether those systems spoil or admit
+that critical behaviour. Among many fascinating
 \cite{jose_renormalization_1977, kankaala_theory_1993,
 ala-nissila_numerical_1994, dierker_consequences_1986, selinger_theory_1988}
-results that emerge from systems with one or more of these fields applied,
-it is predicted that $h_4$ is relevant while $h_6$ is not at some
-sufficiently high temperatures below the Kosterlitz--Thouless point
+results that emerge from systems with one or more of these fields applied, it
+is predicted that $h_4$ is relevant while $h_6$ is not at some sufficiently
+high temperatures below the Kosterlitz--Thouless point
 \cite{jose_renormalization_1977}.
 
 \begin{figure}