\input   mac.tex
\input   moremac.tex
\input peqns.tex
% The spaces inserted after equations and before punctuation marks 
% don't work!

{\bf Abstract}
\medskip
\noindent
We have measured the dynamical critical exponents for the
Swendsen-Wang and the Wolff cluster update algorithms, as well as a 
number of variants of these algorithms, for the $q=2$ and $q=3$ Potts
models in two dimensions. We find that although the autocorrelation
lengths differ considerably between algorithms, the critical exponents
are the same.

\bigskip
\bigskip
{\bf 1. Introduction}
\medskip
\noindent
Cluster update algorithms can greatly reduce critical slowing down in 
computer simulations of spin models. The original idea of Swendsen
and Wang 
\ref\SW{R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. {\bf 58}, 86 (1987)}
(SW) was based on the work of Fortuin and Kasteleyn 
\ref\FK{P.W. Kasteleyn and C.M. Fortuin, J. Phys. Soc. Jpn. Suppl. 
{\bf 26s}, 11 (1969); 
C.M. Fortuin and P.W. Kasteleyn, Physica {\bf 57}, 536 (1972)}
, which relates the Potts model 
\ref\Potts{R.B. Potts, Proc. Camb. Phil. Soc. {\bf 48}, 106 (1952);
F.Y. Wu, Rev. Mod. Phys. {\bf 54}, 235 (1982)}
to a random bond percolation model
\ref\Perc{D. Stauffer, Phys. Rep. {\bf 54}, 1 (1978); 
J.W. Essam, Rep. Prog. Phys. {\bf 43}, 830 (1980)}
. The $q$-state Potts model consists of a lattice of spins $\sigma_{i}$ 
which can take $q$ different values, and whose Hamiltonian is
$$ 
H = K \sum_{<i,j>} \delta_{\sigma_{i}\sigma_{j}} \space .
\eqn\Eqham $$
By introducing bonds between neighboring spins, with probability
$$ 
p(\sigma_i,\sigma_j) = \delta_{\sigma_i\sigma_j} (1 - e^{-K}) \space ,
\eqn\Eqbond $$
clusters of bonded spins are created. The SW algorithm consists of generating
all such clusters, then choosing a random spin value for each cluster and 
assigning it to all the sites in that cluster. In the algorithm of Wolff  
\ref\Wolff{U. Wolff, Phys. Rev. Lett. {\bf 62}, 361 (1989)}
, a site is chosen at random and a single cluster constructed around it. All 
the spins in this cluster are then flipped, i.e. changed to a single new 
Potts state (different from the current state of the cluster's spins) 
which is chosen at random from the $q-1$ possible values.
Note that this is different from the SW update, which may leave the spins 
unchanged in some clusters. These algorithms, which update whole clusters 
of spins at a time, result in a dramatic decrease in critical slowing down 
over local algorithms which update only one spin at a time.

As shown by Wolff 
\ref\WO{U. Wolff, ``Asymptotic Freedom and Mass Generation in the O(3)
Nonlinear Sigma Model'', DESY preprint 89-021, (February 1989)}
, the SW and Wolff algorithms are both special cases of the following general 
cluster update algorithm:

\medskip
\noindent
Step 1. The lattice is decomposed into $N_c$ clusters $c_i$, using the bond
probabilities of equation \Eqbond.

\noindent
Step 2. A subset $C$ of the $c_i$ is chosen using some specified probability 
distribution, and all of the clusters of spins in $C$ are flipped
independently.

\medskip
\noindent
The SW algorithm corresponds to including each cluster in $C$ with probability
$(q-1)/q$, while the Wolff algorithm amounts to picking a single cluster with 
probability proportional to its size. Hence the only differences between the 
Wolff and SW algorithms are that the SW algorithm flips more clusters per 
iteration, while the Wolff algorithm flips relatively larger clusters. In order
to investigate these differences further, we also studied some other methods of
choosing the clusters to be flipped, including the following :
\medskip
\item{(1)}   flip all the clusters (hereafter denoted by AC),
\item{(2)}  flip one of the clusters, chosen at random (denoted by RC), and 
\item{(3)} flip the largest cluster (denoted by LC). 
\medskip
\noindent
Note that flipping all of the clusters will not work in the case of the Ising 
model ($q=2$), since it is equivalent to a global spin flip, and is thus highly 
non-ergodic! However for $q>2$ it is a perfectly acceptable algorithm.

\bigskip
{\bf 2. Measurement of autocorrelations}
\medskip
\noindent
We have made a high statistics study of the above cluster algorithms for the
$q=2$ and $q=3$ Potts models in two dimensions. We measured the 
autocorrelations in the energy, magnetization and susceptibility, and also in 
the improved estimators of these quantities
\Wolff 
\ref\Sweeney{M. Sweeny, Phys. Rev. B {\bf 27}, 4445 (1983)}.
%
%The estimator for the magnetization squared is $<|c|_{SW}^2>$, which for the 
%Wolff algorithm is just the average cluster size $<|c|_{W}>$. The estimator 
%for the energy is
%$$
%E = \sum_{<i,j>} \Theta(i,j) + (1 - \Theta(i,j))/q
%\eqn\Eqimenergy $$ 
%where $\Theta(i,j)$ is ......
%** Sweeney and Wolff argument for why this is **
%
The normalized autocorrelation function for an operator A is
$$ 
\rho(t) = {<A(0)A(t)> - <A(0)>^2 \over <A(0)A(0)> - <A(0)>^2} \space . 
\eqn\Eqautocorr $$
At large t this should go as $e^{-t/\tau_{exp}}$, where $\tau_{exp}$ is the 
exponential autocorrelation length. The quantity which is actually relevant 
to the error in A is the integrated autocorrelation length 
\ref\Sokal{N. Madras and A.D. Sokal, J. Stat. Phys. {\bf 50}, 109 (1988)}
$$ 
\tau_{int} = {1 \over 2} + \sum_{t=1}^{\infty} \rho(t) \space . 
\eqn\Eqintauto $$
%If $\rho(t)$ is exactly exponential $\tau_{int} = \tau_{exp}$ to a good 
%approximation. 
Both the integrated and exponential autocorrelation lengths have an 
associated exponent $z$ governing critical slowing down given by
$\tau(L) \sim L^{z}$, where $L$x$L$ is the size of the lattice and $\tau(L)$
is measured at the critical point.

Autocorrelations are traditionally measured between each update of the
entire lattice, so for the single cluster updates, where only a fraction of
the lattice sites are updated at each iteration, the autocorrelation length
needs to be scaled by a factor of the average cluster size divided by the 
lattice size, so that
$$  
\tau = \tau' <|c|> L^{-2} \space ,
\eqn\Eqscale $$
where $\tau'$ is the unscaled autocorrelation length.
For the Wolff algorithm, the average cluster size scales like the 
susceptibility, hence $<|c|_{W}> \sim L^{\gamma/\nu}$. The average cluster 
size for SW seems to be approximately independent of $L$, thus
$<|c|_{RC}> \sim L^0$. 
As in the percolation model \Perc,
the maximum cluster size (normalized to the size of the lattice) scales 
like the magnetization, so that $<|c|_{LC}> \sim L^{2-\beta/\nu}$. 
From equation \Eqscale\ and the definition of $z$ we have that
$$ \eqalign{
z_{W} &= z_{W}' - (2 - \gamma / \nu) \cr
z_{RC} &= z_{RC}' - 2 \cr
z_{LC} &= z_{LC}' - \beta/\nu \cr
}\eqn\Eqzscale $$
where $z'$ is the unscaled exponent, $\tau'(L) \sim L^{z'}$. In our 
calculations the scaled value of $z$ is actually obtained by scaling the 
autocorrelations by the average cluster size, but we have checked that
using equation \Eqzscale\ instead gives consistent answers for LC and
Wolff. For RC, however, there are large finite size effects on the smaller 
lattices and the average cluster size does not scale well.
Note that the computational requirement for the Wolff algorithm is 
proportional to the average size of the clusters which are flipped \Wolff, 
but for the other single cluster algorithms (RC and LC) it is still
necessary to form all the clusters in order to choose the one to flip.
Hence the computational requirement goes like $L^2$ (the same as for SW), 
and in that sense 
the above scaling is unwarranted. However it may be possible to invent
a clever algorithm which constructs single clusters with the required
probability distribution, and in any case the scaling is necessary to
allow a sensible comparison with the Wolff algorithm.

The autocorrelations were measured at the infinite volume critical point
(which is known
exactly for 2-$d$ Potts models) over runs of typically $10^5$ to $10^6$ 
sweeps ($10^6$ to $10^7$ sweeps for the RC algorithm),
and 10 or 20 such runs were done for each of the different lattice 
sizes, from $8^2$ up to $128^2$. The results 
from the multiple runs were averaged to give the overall autocorrelation
function $\rho(t)$. 
For the SW algorithm, and for the Wolff algorithm in the case $q=2$, we 
also did a single run of $10^6$ sweeps on a $256^2$ lattice.
The autocorrelation length $\tau_{exp}$ was extracted by a 
fit to $log$ $\rho(t)$ in a region where it was linear, typically 
somewhere between $\tau_{exp}$ and $3\tau_{exp}$.
%
%** Say something about Wolff not giving good fits ?? **
%
$\tau_{int}$ was calculated by summing $\rho(t)$ for $t$ less than some
cutoff value, 
and using the exponential fit to $\rho(t)$ for the remaining part of 
the sum, which is then just a geometric series and so can be summed 
analytically. The cutoff was taken to be the largest point used in the 
exponential fit, which was usually between $2\tau_{exp}$ and $3\tau_{exp}$. 
This method is similar to that used by Wolff
\ref\WIsing{U. Wolff, Phys. Lett. {\bf B228}, 379 (1989)}.
The above analysis was also done for each of the 10 or 20 multiple runs, which 
provided an estimate of the error in the autocorrelation lengths.

\bigskip
{\bf 3. Results}
\medskip
\noindent
We found that the energy, magnetization and susceptibility all have very
similar values of $\tau_{exp}$, however the global quantities (especially
the magnetization) decorrelate faster than the energy, giving smaller values
of $\tau_{int}$, and smaller values of $z_{int}$ as well.
For the SW algorithm the improved estimators gave approximately the same
autocorrelations as the standard values of the measured quantities. 
For the Wolff algorithm, the improved estimator of the (unsubtracted)
susceptibility is just the average cluster size. 
This quantity seems to have approximately the same 
$z_{exp}$ as the standard estimators, although $z_{int}$ is very much 
reduced. In fact for $q=2$ $\tau_{int}$ actually decreases slightly after 
the scaling of equation \Eqscale.
Unfortunately this does not mean that the improved 
estimator gives a greatly reduced error in the susceptibility, in fact at 
criticality the error is greater than that of the standard estimator due to 
the very large fluctuations in size of the clusters which are grown.

Results for $\tau_{int}$ of the energy for the different algorithms are
shown in figs. 1 and 2 for $q=2$ and $q=3$ respectively, along with 
straight lines which represent chi-squared fits to a 
power law. The values of $z$ from fits to all the data are given in Table 1.
We have not given any values for the RC algorithm, since the data is not 
good enough to reliably extract an exponent.
The scaled exponents for all the other algorithms are approximately the same, 
and in particular it appears that $z_{W} = z_{SW}$ 
for Potts models in two dimensions. The autocorrelation lengths differ 
substantially between algorithms, however. The Wolff algorithm has values
which are less than those of SW by a factor of about 0.45 for $q=2$ and 0.67 
for $q=3$.
It is interesting to note that these ratios are practically the same as the 
probabilities for flipping a cluster in the SW algorithm,
so that if the SW algorithm was scaled (like the single cluster algorithms) by 
the proportion of spins which are actually {\it changed} in the update, then 
the Wolff and SW algorithms would give nearly identical results.
It therefore appears that the main advantage of a single 
cluster algorithm such as Wolff is that the cluster is always flipped.
Notice however that flipping all the SW clusters (the AC algorithm) is
actually worse than the standard SW algorithm. 
Another proposed reason for the superiority of the Wolff algorithm over 
SW is that the average cluster size is larger, and this is supported by the
fact that the LC algorithm, where the largest cluster is flipped,
does slightly better than Wolff. It is therefore quite surprising to see that
the RC algorithm, where the average cluster size is approximately the same as
for SW, does even better than Wolff or LC. This suggests that (at least in
two dimensions) it is the single cluster nature of these updates, rather than 
the size of the clusters, which is most important. 

We have also done preliminary investigations of two other possible update
algorithms, with a single run of $10^5$ to $10^6$ sweeps for each lattice 
size. The first algorithm is just SW except that the largest cluster 
is always flipped. It gave autocorrelation lengths which agreed with those 
of SW within errors. In the second algorithm, all ``large'' clusters (those 
greater than a certain cutoff size) are flipped. In a study of the 
4-$d$ mean field Ising model, this algorithm was found to give the same 
value of $z$ as SW, but the values of $\tau$ were halved
\ref\MFI{T.S. Ray, P. Tamayo and W. Klein, Phys. Rev. A {\bf 39}, 5949 (1989)}.
For the 2-$d$ Potts model we chose the cutoff
to be $L^{3/2}$, and found that the autocorrelation lengths for this 
algorithm are decreased only slightly, giving values approximately midway
between those of SW and LC, with the same value of $z$.

While this work was in progress a number of other results have appeared 
concerning the dynamical exponents of cluster algorithms for Potts models. 
Wolff \WIsing\
has compared his algorithm with SW for the 2-$d$ and 
3-$d$ Ising models. At larger lattice sizes we find slightly higher values 
of $\tau$ than he quotes, along with correspondingly higher values of $z$,
however we still find a value of $z_{SW}$ which is somewhat less than that 
obtained by Swendsen and Wang \SW.
A study of the Wolff algorithm for the Ising model has also been carried
out by Tamayo et al.
\ref\Brower{P. Tamayo, R.C. Brower and W. Klein, ``Single Cluster Monte Carlo 
Dynamics for the Ising Model'', Boston University preprint, (July 1989)}
, who obtain similar results. In particular they show that $z_{W} < z_{SW}$ in 
three and four dimensions. Their paper also contains some remarks on 
the RC algorithm, and suggests that $z_{RC} < z_{W}$ for the 2-$d$ model.
From our results it 
appears that this is certainly possible, however the trend at larger lattice 
sizes seems to imply that the exponent is probably the same as for the 
other algorithms.
Much better data would be required to draw a firm conclusion.
Our value for $z_{SW}$ in the $q = 3$ case agrees with that of Swendsen 
and Wang \SW, and the more recent result of Li and Sokal
\ref\Li{X.-J. Li and A.D. Sokal, Phys. Rev. Lett. {\bf 63}, 827 (1989)}.
Kandel et al.
\ref\Kandel{D. Kandel, E. Domany and A. Brandt, Phys. Rev. B {\bf 40}, 
330 (1989)}
obtain a value of 0.4 by measuring the relaxation time from an ordered state 
to equilibrium, claiming to see a crossover to this smaller value of the 
exponent for $L > 32$, however we see no sign of such behaviour.  
%
% It is possible that the different results are due to a cluster 
% dynamics equivalent of the difference between the equilibrium and 
% non-equilibrium exponents for local update algorithms
% \ref\Racz{Z. Racz, Phys. Rev. B {\bf 13} 263 (1976); M.E. Fisher and 
% Z. Racz, Phys. Rev. B {\bf 13} 5039 (1976)}.
%
An interesting recent development is the claim by Burkitt and Heermann
that the autocorrelations in the 2-$d$ Ising model using the SW algorithm 
grow logarithmically rather than as a power law
\ref\Burkitt{D.W. Heermann and A.N. Burkitt, ``System Size Dependence of 
the Autocorrelation Time for the Swendsen-Wang Ising Model'',
Bergische University preprint, (1989)}
, so that $z_{SW} = 0$. 
From our data it is difficult to distinguish between a power law and a
logarithm. 
If we look at the slopes of lines connecting successive points in fig. 1,
we do see a slight trend towards smaller values as $L$ increases, which is
not evident in the $q=3$ data, where a logarithmic increase in $\tau$ can be
definitely excluded. However it is not clear that the $q=2$ data does not
follow a power law, especially since the statistics for 
the largest lattice sizes are not particularly good.
We are currently collecting more data to try to resolve this point.  

\bigskip
{\bf Acknowledgements}
\medskip
\noindent
The simulations were done using a 32 processor Meiko Computing Surface,
a 20 processor Sequent Symmetry, a 20 processor Encore Multimax,
a 96 processor BBN GP1000 Butterfly and a network of SUN workstations.
We would like to thank the Caltech Concurrent Supercomputer Facility
and the Advanced Computing Research Facility at Argonne National Laboratory
for the use of these machines.
This work was sponsored in part by DOE grants
DE-FG03-85ER25009 and DE-AC03-81ER40050.

\listrefs\vfill
\end
