% Created by Jaron Kent-Dobias on Thu Apr 20 12:50:56 EDT 2017. % Copyright (c) 2017 Jaron Kent-Dobias. All rights reserved. \documentclass[aps,pre,reprint]{revtex4-1} \usepackage{amsmath,amssymb,latexsym,mathtools,algorithm,listings,hyperref} \hypersetup{breaklinks=true} % uncomment to label only equations that are referenced in the text %\mathtoolsset{showonlyrefs=true} % I want equation labels but don't want to type out `equation' \def\[{\begin{equation}} \def\]{\end{equation}} % additional math \def\tr{\mathop{\mathrm{Tr}}\nolimits} % trace \def\sgn{\mathop{\mathrm{sgn}}\nolimits} % sign function \def\dd{d} % differential \def\Z{\mathbb Z} % integers \def\R{\mathbb R} % reals % physics conventions \def\c{\mathrm c} % critical label as in T_c \def\H{\mathcal H} % Hamiltonian \def\J{Z} % site-site coupling \def\B{B} % external field \newcommand\set[1]{\mathbf #1} % collection of N spins \newcommand\avg[1]{\langle#1\rangle} % ensemble average % dimensions \def\dim{d} \def\twodee{\textsc{2\dim} } \def\threedee{\textsc{3\dim} } \def\fourdee{\textsc{4\dim} } % fancy partial derivative \newcommand\pd[3][]{ \ifthenelse{\isempty{#1}} {\def\tmp{}} {\def\tmp{^#1}} \frac{\partial\tmp#2}{\partial#3\tmp} } % C++, but pretty \def\CXX{C\texttt{++}} \begin{document} \title{Cluster representations and the Wolff algorithm in arbitrary external fields} \author{Jaron Kent-Dobias} \author{James P.~Sethna} \affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA} \date\today \begin{abstract} We introduce a natural way to extend celebrated spin-cluster Monte Carlo algorithms for fast thermal lattice simulations at criticality, like Wolff, to systems in arbitrary fields, be they linear magnetic vector fields or nonlinear anisotropic ones. 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 which lives in the space of symmetry transformations. The ordinary cluster-building process can then be run on the new representation. We show that this extension preserves the scaling of accelerated dynamics in the absence of a field for Ising, Potts, and $\mathrm O(n)$ models and demonstrate the method's use in modelling the presence of novel nonlinear fields. We also provide a \CXX\ library for the method's convenient implementation for arbitrary models. \end{abstract} \maketitle Lattice models are important in the study of statistical physics and phase transitions. Rarely exactly solvable, they are typically studied by approximate and numerical methods. Monte Carlo techniques are a common way of doing this, approximating thermodynamic quantities by sampling the distribution of 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. Celebrated cluster algorithms largely addressed this in the absence of symmetry-breaking fields by using nonlocal updates \cite{janke_nonlocal_1998} whose clusters undergo a percolation transition at the critical point of the system \cite{coniglio_clusters_1980}. These result in relatively small dynamic exponents for many spin systems \cite{wolff_comparison_1989, du_dynamic_2006, liu_dynamic_2014, wang_cluster_1990}, including the Ising, $\mathrm O(n)$ \cite{wolff_collective_1989}, and Potts \cite{swendsen_nonuniversal_1987, baillie_comparison_1991} models. These algorithms rely on the natural 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 returns global rotation invariance to spin Hamiltonians at the cost of an extra degree of freedom, allowing the method to be used in a subcategory of interesting fields \cite{alexandrowicz_swendsen-wang_1989, wang_clusters_1989, ray_metastability_1990}. Static fields have also been applied by including a separate metropolis or heat bath update step after cluster formation \cite{destri_swendsen-wang_1992, lauwers_critical_1989, 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 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 several models suggests that the `ghost' approach is a natural one, e.g., that it extends the celebrated scaling of dynamics in these algorithms at zero field to various non-symmetric perturbations. We also show, by a redefinition of the spin--spin coupling in a generic class of spin systems, \emph{arbitrary} external fields can be treated using cluster methods. Rather than the introduction of a `ghost spin,\!' our representation relies on introducing a `ghost transformation,\!' 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. We provide an open-source implementation of this method in the form of a \CXX\ library, available at \url{https://git.kent-dobias.com/wolff/} \cite{kent-dobias_wolff_2018}. Use of this library will be described briefly within, but extensive documentation is also available at \url{https://doc.kent-dobias.com/wolff/}. \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 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 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 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 \[ \H(\set s)=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!\J(s_i,s_j)-\sum_{i\in V}\B(s_i), \] where $\J:X\times X\to\R$ couples adjacent spins and $B:X\to\R$ is an external field. $\J$ 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$, $\J(r\cdot s,r\cdot t)=\J(s,t)$. One may also allow $\J$ to also be a function of edge---for modelling random-bond, long-range, or anisotropic interactions---or allow $B$ to be a function of site---for applying arbitrary boundary conditions or modelling random fields. The formal results of this paper (that the algorithm obeys detailed balance and ergodicity) hold equally well for these cases, but we will drop the additional index notation for clarity. Some extensions, like adding strong random fields or bonds, ultimately prove inefficient \cite{rieger_monte_1995, redner_graphical_1998}. \begin{table*}[htpb] \begin{tabular}{l||ccccc} & Spins ($X$) & Symmetry ($R$) & Action ($g\cdot s$) & Coupling ($\J(s,t)$) & Common Field ($B(s)$) \\ \hline\hline Ising & $\{-1,1\}$ & $\Z/2\Z$ & $0\cdot s\mapsto s$, $1\cdot s\mapsto -s$ & $st$ & $Hs$ \\ $\mathrm O(n)$ & $S^{n-1}$ & $\mathrm O(n)$ & $M\cdot s\mapsto Ms$ & $s^{\mathrm T}t$ & $H^{\mathrm T}s$\\ Potts & $\{1,\ldots,q\}$ & $\mathrm S_n$ & $(i_1,\ldots,i_q)\cdot s=i_s$ & $\delta(s,t)$ & $\sum_mH_m\delta(m,s)$\\ Clock & $\Z/q\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)$\\ \textsc{Dgm} & $\Z$ & $D_{\mathrm{inf}}$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $(s-t)^2$ & $Hs^2$\\ \end{tabular} \caption{Several examples of spin systems and the symmetry groups that act on them. Common choices for the spin--spin coupling in these systems and their external fields are also given. Other fields are possible, of course: for instance, some are interested in modulated fields $H\cos(2\pi k\theta(s))$ for integer $k$ and $\theta(s)$ giving the angle of $s$ to some axis applied to the $\mathrm O(2)$ model \cite{jose_renormalization_1977}. All models listed here have example implementations in the provided \CXX\ library \cite{kent-dobias_wolff_2018}. } \label{table:models} \end{table*} Implementation of a model in the provided library is as simple as defining a class that represents an element of the state space $X$, with default constructor (and destructor, if necessary), and a class that represents an element of the group $R$, with default constructor and member functions that define the action and inverse action of the class on both states and group elements. Specific details may be found at \url{https://doc.kent-dobias.com/wolff/models.html}. The goal of statistical mechanics is to compute expectation values of observables $A:X^N\to\R$. Assuming the ergodic hypothesis holds (for systems with broken-symmetry states, it does not), the expected value $\avg A$ of an observable $A$ is its average over every state $\set s$ in the configuration space $X^N$ weighted by the Boltzmann probability of that state appearing, or \[ \avg A =\frac{\int_{X^N}A(\set s)e^{-\beta\H(\set s)}\,\dd\mu(\set s)} {\int_{X^N}e^{-\beta\H(\set s)}\,\dd\mu(\set s)}, \] where for $Y_1\times\cdots\times Y_N=Y\subseteq X^N$ the product measure $\mu(Y)=\mu(Y_1)\cdots\mu(Y_N)$ is the simple extension of the measure on $X$ to a measure on $X^N$. These values are estimated using Monte Carlo techniques by constructing a finite sequence of states $\{\set{s_1},\ldots,\set{s_M}\}$ such that \[ \avg A\simeq\frac1M\sum_{i=1}^MA(\set{s_i}). \] Sufficient conditions for this average to converge to $\avg A$ as $M\to\infty$ are that the process that selects $\set{s_{i+1}}$ given the previous states be Markovian (only depends on $\set{s_i}$), ergodic (any state can be accessed), and obey detailed balance (the ratio of probabilities that $\set{s'}$ follows $\set s$ and vice versa is equal to the ratio of weights for $\set s$ and $\set{s'}$ in the ensemble). Measurements of observables during Monte Carlo in the provided library are made by the use of hooks, which are member functions of a measurement class that are run at designated points during the algorithm's execution and are provided arbitrary information about the internal state of all relevant objects. A detailed description of these hooks can be found at \url{https://doc.kent-dobias.com/wolff/measurement.html}. While any of several related cluster algorithms can be described for this system, we will focus on the Wolff algorithm \cite{wolff_collective_1989}. In the absence of an external field, e.g., $B(s)=0$, the Wolff algorithm proceeds in the following way.\\ \begin{algorithm}[H] \begin{enumerate} \item Pick a random site $m_0$ and add it to the stack. \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 set of spin states is infinite. \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(\J(r\cdot s_m,s_j)-\J(s_m,s_j))}\}. \label{eq:bond_probability} \] \item Take $s_m\mapsto r\cdot s_m$. \end{enumerate} \end{enumerate} \end{enumerate} \caption{Wolff} \label{alg:wolff} \end{algorithm} When the stack is exhausted, a cluster of connected spins will have been 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 $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 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} \frac{P(\set s\to\set{s'})}{P(\set{s'}\to\set s)} &=\frac{f(r\mid m_0,\set s)}{f(r^{-1}\mid m_0,\set s')} \prod_{\{i,j\}\in }\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\partial C}e^{\beta(\J(s_i',s_j')-\J(s_i,s_j))} =\frac{e^{-\beta\H(\set{s'})}}{e^{-\beta\H(\set s)}}, \end{aligned} \] \end{widetext} whence detailed balance is also satisfied, using $r=r^{-1}$ and $Z(r\cdot 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 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, caracciolo_wolff-type_1993}. This algorithm can be run on a system using the provided library. To construct a system, you must provide a graph representing the lattice, a temperature, the spin coupling function $Z$, and the field coupling function $B$. Once constructed, cluster flips as described in Alg.~\ref{alg:wolff} can be performed by directly providing seed sites $m_0$ and transformations $r$, or many in sequence by providing a function that generates random (appropriately distributed to preserve detailed balance) transformations $r$. The construction and use of Wolff systems is described at \url{https://doc.kent-dobias.com/wolff/system.html}. \section{Adding the field} This algorithm relies on the fact that the coupling $\J$ depends only on relative orientation of the spins---global reorientations do not affect the Hamiltonian. The external field $B$ breaks this symmetry. Fortunately it can be restored. Define a new graph $\tilde G=(\tilde V,\tilde E)$, where $\tilde V=\{0,1,\ldots,N\}$ adds the new `ghost' site $0$ which is connected by \[ \tilde E=E\cup\big\{\{0,i\}\mid i\in V\big\} \] to all other sites. Instead of assigning the ghost site a spin whose value comes from $X$, we assign it values in the symmetry group $s_0\in R$, so that the configuration space of the new model is $R\times X^N$. We introduce the Hamiltonian $\tilde\H:R\times X^N\to\R$ defined by \[ \begin{aligned} \tilde\H(s_0,\set s) &=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!\J(s_i,s_j) -\sum_{i\in V}B(s_0^{-1}\cdot s_i)\\ &=-\!\!\!\!\sum_{\{i,j\}\in\tilde E}\!\!\!\!\tilde\J(s_i,s_j), \end{aligned} \] where the new coupling $\tilde\J:(R\cup X)\times(R\cup X)\to\R$ is defined for $s,t\in R\cup X$ by \[ \tilde\J(s,t) = \begin{cases} \J(s,t) & \text{if $s,t\in X$} \\ B(s^{-1}\cdot t) & \text{if $s\in R$} \\ B(t^{-1}\cdot s) & \text{if $t\in R$}. \end{cases} \label{eq:new.z} \] The modified coupling is invariant under the action of group elements: for any $r,s_0\in R$ and $s\in X$, \[ \begin{aligned} \tilde\J(rs_0,r\cdot s) &=B((rs_0)^{-1}\cdot (r\cdot s))\\ &=B(s_0^{-1}\cdot s) =\tilde\J(s_0,s) \end{aligned} \] 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 old system, whose properties we actually want to measure? If $A:X^N\to\R$ is an observable of the original system, we construct an observable $\tilde A:R\times X^N\to\R$ of the new system defined by \[ \tilde A(s_0,\set s)=A(s_0^{-1}\cdot\set s) \] whose expectation value in the new system equals that of the original observable in the old system. First, note that $\tilde\H(1,\set s)=\H(\set s)$. Since the Hamiltonian is invariant under global rotations, it follows that for any $g\in R$, $\tilde\H(g,g\cdot\set s)=\H(\set s)$. Using the invariance properties of the measure on $X$ and introducing a measure $\rho$ on $R$, it follows that \[ \begin{aligned} \avg{\tilde A} &=\frac{ \int_R\int_{X^N}\tilde A(s_0,\set s)e^{-\beta\tilde\H(s_0,\set s)}\,\dd\mu(\set s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\set s)}\,\dd\mu(\set s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(s_0^{-1}\cdot\set s)e^{-\beta\tilde\H(s_0,\set s)}\,\dd\mu(\set s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\set s)}\,\dd\mu(\set s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(\set{s'})e^{-\beta\tilde\H(s_0,s_0\cdot\set{s'})}\dd\mu(s_0\cdot\set{s'})\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,s_0\cdot\set{s'})}\dd\mu(s_0\cdot\set{s'})\,\dd\rho(s_0) }\\ &=\frac{ \int_R\dd\rho(s_0)}{ \int_R\dd\rho(s_0)}\frac{\int_{X^N}A(\set{s'})e^{-\beta\H(\set{s'})}\dd\mu(\set{s'}) }{\int_{X^N}e^{-\beta\H(\set{s'})}\dd\mu(\set{s'}) } =\avg A. \end{aligned} \] Using this equivalence, spin systems in a field may be treated in the following way. \begin{enumerate} \item Add a site to your lattice adjacent to every other site. \item Initialize a `spin' at that site whose value is a representation of a member of the symmetry group of your ordinary spins. \item Carry out the ordinary Wolff cluster-flip procedure on this new lattice, substituting $\tilde\J$ as defined in \eqref{eq:new.z} for $\J$. \end{enumerate} Ensemble averages of observables $A$ can then be estimated by sampling the value of $\tilde A$ on the new system. In contrast with the simpler ghost spin representation, this form of the Hamiltonian might be considered the `ghost transformation' representation. 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, Potts, and (with clusters weighted by the components of their spins along the reflection direction \cite{hasenbusch_improved_1990}) $\mathrm O(n)$ models, 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 and the squared magnetization continues to hold (see \eqref{eq:cluster-size-scaling} and Fig.~\ref{fig:cluster_scaling} below). \section{Examples} \label{sec:examples} Several specific examples from Table~\ref{table:models} are described in the following. \subsection{The Ising model} In the Ising model spins are drawn from the set $\{1,-1\}$. Its symmetry group is $C_2$, the cyclic group on two elements, which can be conveniently represented by a multiplicative group with elements $\{1,-1\}$, exactly the same as the spins themselves. The only nontrivial element is of order two, and is selected every time in the algorithm. Since the symmetry group and the spins are described by the same elements, performing the algorithm on the Ising model in a field is fully described by just using the `ghost spin' representation. This algorithm or algorithms based on the same decomposition of the Hamiltonian have been applied by several researchers \cite{alexandrowicz_swendsen-wang_1989, wang_clusters_1989, ray_metastability_1990}. The algorithm has been implemented by one of the authors in an existing interactive Ising simulator at \url{https://mattbierbaum.github.io/ising.js} \cite{bierbaum_ising.js_2016}. \subsection{The $\mathrm O(n)$ models} \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 $\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, taking advantage of the freedom to choose the function $f$ in Alg.~\ref{alg:wolff}. 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, 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 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 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 comparable to the lattice spacing and is efficiently simulated 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} \caption{ The scaled autocorrelation time of the energy $\H$ for the Wolff algorithm on a $32\times32\times32$ \textsc{xy} model at its critical temperature as a function of applied vector field magnitude $|H|$. Red points correspond to reflections sampled uniformly, while the green points represent reflections sampled as described in section \ref{sec:examples:on}. } \label{fig:generator_times} \end{figure} \subsection{The Potts model} In the $q$-state Potts model spins are described by elements of $\{1,\ldots,q\}$. Its symmetry group is the symmetric group $\mathrm S_n$ of permutations of its elements. The element $(i_1,\ldots,i_q)$ takes the spin $s$ to $i_s$. There are potentially many elements of order two, but the two-element swaps alone are sufficient to both generate the group and act transitively on $\{1,\ldots,q\}$, providing ergodicity. \subsection{Clock models} In the $q$-state clock model spins are described by elements of $\Z/q\Z$, the set of integers modulo $q$. Its symmetry group is the dihedral group $D_q=\{r_0,\ldots,r_{q-1},s_0,\ldots,s_{q-1}\}$, the group of symmetries of a regular $q$-gon. The element $r_n$ represents a rotation by $2\pi n/q$, and the element $s_n$ represents a reflection composed with the rotation $r_n$. The group acts on spins by permutation: $r_n\cdot m={n+m}\pmod q$ and $s_n\cdot m={-(n+m)}\pmod q$. This is the natural action of the group on the vertices of a regular polygon that have been numbered $0$ through $q-1$. The elements of $D_q$ of order 2 are all reflections and $r_{q/2}$ if $q$ is even, though the former can generate the latter. While reflections do not necessarily generate the entire group, their action on $\Z/q\Z$ is transitive and therefore the algorithm is ergodic. \subsection{Roughening models} Though not often thought of as a spin model, roughening of surfaces can be described in this framework. Spins are described by integers $\Z$ and their symmetry group is the infinite dihedral group $D_\infty=\{r_i,s_i\mid i\in\Z\}$, whose action on the spin $j\in\Z$ is given by $r_i\cdot j=i+j$ and $s_i\cdot j=-i-j$. The elements of order two are reflections $s_i$, whose action on $\Z$ is transitive. The coupling can be any function of the absolute difference $|i-j|$. Because 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 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} are modified by adding a constant $0