\input /net/jack/paulc/texlib/mac.tex
\input /net/jack/paulc/texlib/moremac.tex
% \input /net/jack/paulc/texlib/fonts.tex
\input /net/jack/paulc/texlib/peqns.tex

\font\big=cmr10 scaled\magstep1
\font\bigger=cmr10 scaled\magstep2

\line{SCCS-???\hfill}
\line{CRPC-TR91???\hfill}
\line{COLO-HEP-???\hfill}
\bigskip
\bigskip
\bigskip
\centerline{\bf \bigger{DYNAMICAL CRITICAL EXPONENTS FOR}}
\smallskip
\centerline{\bf \bigger{ISING MODEL CLUSTER ALGORITHMS}}
\bigskip
\bigskip
\bigskip
\centerline{\big Clive F. BAILLIE} 
\medskip
\centerline{\it Physics Department, University of Colorado, Boulder CO 80309, U.S.A.}
\bigskip
\centerline{\big Paul D. CODDINGTON and Leping HAN}
\medskip
\centerline{\it Physics Department, Syracuse University, Syracuse NY 13244, U.S.A.}
\bigskip
\bigskip
\centerline{September, 1991}
\bigskip
\bigskip
\bigskip
\centerline{\bf Abstract}
\medskip
\noindent
We have measured the dynamical critical exponents for the
Swendsen-Wang and the Wolff cluster update algorithms for the Ising model 
in 2, 3, and 4 dimensions.
The data for the Wolff algorithm suggest that 
$tau_{int,E} = constant \times C_{H}$,
so that $z_{int,E} = \alpha / \nu$.
We also compare our data with results for models in the same universality 
class as the Ising models, to investigate dynamical universality for 
cluster algorithms.
%
%\bigskip
%\bigskip
%\noindent
%PACS numbers 05.50.+q, 11.15.Ha, 64.60.Ht.
%
\vfill
\eject

\bigskip
{\bf 1. \quad Introduction}
\medskip

The Monte Carlo cluster update algorithms of Swendsen and Wang (SW)
\ref\SW{R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. {\bf 58} (1987) 86}
and Wolff
\ref\Wolff{U. Wolff, Phys. Rev. Lett. {\bf 62} (1989) 361}
can greatly reduce critical slowing down in computer simulations of spin 
models. 
Since they were introduced, a large amount of work has gone into
studying these algorithms, 
and in particular, in measuring their dynamical
critical exponents for various spin models, in an attempt to understand 
how these algorithms work, 
and how they can be generalized to more complicated models
(for reviews of cluster algorithms, see refs.
\ref\Sokalrev{A.D. Sokal, in: Computer Simulation Studies in Condensed Matter
Physics:  Recent Developments, eds. D.P. Landau, K.K. Mon
and H.-B. Sch\"uttler (Springer-Verlag, Berlin-Heidelberg, 1988);
A.D. Sokal, Proc. of the International Conference on Lattice Field Theory,
Tallahassee, (Oct. 1990), Nucl. Phys. B (Proc. Suppl.). {\bf 20} (1991) 55},
\ref\reviews{U. Wolff, Proc. of the International Conference on Lattice Field
Theory, Capri, (Sept. 1989), Nucl. Phys. B (Proc. Suppl.) {\bf 17} (1990) 93;
J.-S. Wang and R.H. Swendsen, Physica A {\bf 167} (1990) 565;
C.F. Baillie, Int. J. of Mod. Phys. C, {\bf 1} (1990) 91}
).
%
% More waffle?
%

Experience has shown that it is very difficult to calculate
dynamical exponents, and that very good statistics and very large lattices
are required to get good results, for both cluster algorithms and local
algorithms
\ref\zMetrop{
% R.B. Pearson, J.L. Richardson and D. Toussaint, Phys. Rev. B {\bf 31}
% (1985) 4472;
% N. Ito, M. Taiji and M. Suzuki, J. Phys. Soc. Japan {\bf 56} (1987) 4218;
S. Tang and D.P. Landau, Phys. Rev. B {\bf 36} (1987) 567;
S. Wansleben and D.P. Landau, Phys. Rev. B {\bf 43} (1991) 6006}.
As an example of the difficulties involved, consider the case of the 
two dimensional Ising model
\ref\Ising{B.M. McCoy and T.T. Wu, The Two-Dimensional Ising
Model, (Harvard University Press, Cambridge, Mass., 1973);
F.Y. Wu, Rev. Mod. Phys. {\bf 54} (1982) 235}.
Swendsen and Wang estimated $z$ for their algorithm to be 0.35(1) \SW.
Tamayo {\it et al.} 
\ref\Tamayo{P. Tamayo, R.C. Brower and W. Klein,
J. Stat. Phys. {\bf 58} (1990) 1083}
obtained 0.33(5) for the Wolff algorithm.
More accurate measurements by Wolff 
\ref\WIsing{U. Wolff, Phys. Lett. {\bf B228} (1989) 379}
gave 0.27(2) and 0.26(2) for the dynamical exponent of the energy for the
SW and Wolff algorithms respectively, and 0.20(2) and 0.13(2) for the
exponent of the susceptibility.
Subsequent work on larger lattices by Burkitt and Heermann 
\ref\Burkitt{D.W. Heermann and A.N. Burkitt, Physica A {\bf 162} (1990) 210}
claimed that the autocorrelation times at the critical
point grew logarithmically with lattice size, so that $z = 0$.
Recent results suggest that it is very difficult to distinguish between
a logarithm and a small power, even with long runs on lattices up to $512^2$,
especially if there are non-negligible corrections to scaling
\ref\us{C.F. Baillie and P.D. Coddington, Phys. Rev. B {\bf 43} (1991) 10617}.
These results are summarized in Table 1.

Measurements for the three dimensional Ising model have proven to be just as
difficult. 
For the SW algorithm, the original measurement by Swendsen and Wang yielded 
a value of 0.75(1) \SW; 
a later calculation by Wolff gave a value of 0.50(3) for both the energy
and susceptibility \WIsing;
while recent work by Wang gave the value of 0.46(4) and 0.339(4) for different
magnetization operators, with the claim that the discrepancy may be due to
corrections to scaling
\ref\Wang{J.-S. Wang, Physica A {\bf 164} (1990) 240}.
For the Wolff algorithm, Tamayo {\it et al.} \Tamayo
obtained 0.44(10);
Wolff used better data to obtain values of 0.28(2) for the energy exponent and 
0.14(2) for the susceptibility exponent \WIsing.
We have subsequently re-examined Wolff's data, and found that 
it also fits well to a logarithm, so that $z = 0$ is also a possibility for 
the Wolff algorithm.

In four dimensions there is only one result known, which is $z = -0.05(15)$
for the Wolff algorithm \Tamayo. 
Simulations have also been done on the mean-field Ising model, which is
equivalent to an infinite dimensional Ising model.
%
% Should it give the same answer as d >= 4 Ising model (or d > 4) ???
%
These give values of -0.10(15) for the Wolff algorithm \Tamayo,
and 0.98(8) for SW,
which agrees with a theoretical argument that $z = 1$ for SW in
mean-field theory
\ref\MF{T.S. Ray, P. Tamayo and W. Klein, Phys. Rev. A {\bf 39} (1989) 5949}.

Due to the discrepancies between the various measurements of the dynamical
critical exponents, we decided to do some simulations with good statistics,
and on fairly large lattices, to try to get reliable values for $z$.
Accurate values for the dynamical exponents are needed in order to be of
help in constructing and checking theories of how cluster algorithms work,
and to compare with conjectured
\ref\Klein{W. Klein, T. Ray and P. Tamayo, Phys. Rev. Lett. {\bf 62} (1989) 163}
and rigorous 
\ref\Li{X.-J. Li and A.D. Sokal, Phys. Rev. Lett. {\bf 63} (1989) 827}
bounds on the exponents.

\bigskip
{\bf 2. \quad Simulations}
\medskip

We have made a high statistics study of the SW and Wolff algorithms for the
Ising model in 2, 3, and 4 dimensions. 
We measured the autocorrelations in the energy 
% $E = - \sum_{<i,j>} \sigma_{i}\sigma_{j}$
and the absolute magnetization
% $|M| = \vert \sum_j \sigma_j \vert$
. The normalized autocorrelation function for an operator $A$ is
$$ 
\rho_A(t) = {<A(0)A(t)> - <A(0)>^2 \over <A(0)A(0)> - <A(0)>^2}. 
\eqn\Eqautocorr
$$
At large $t$ this should decay as $e^{-t/\tau_{exp}}$, where $\tau_{exp}$ 
is the exponential autocorrelation time, which is independent of $A$ 
unless $A$ is orthogonal to the slowest mode of the system
\Sokalrev.
The quantity which is actually relevant 
to the error in $A$ is the integrated autocorrelation time 
\ref\Sokal{N. Madras and A.D. Sokal, J. Stat. Phys. {\bf 50} (1988) 109}
,\Sokalrev
$$ 
\tau_{int,A} = {1 \over 2} + \sum_{t=1}^{\infty} \rho_A(t). 
\eqn\Eqintauto
$$
$\tau_{int,A}$ and $\tau_{exp}$  have associated 
dynamical critical exponents $z_{int,A}$ and $z_{exp}$, which are given by
$\tau(L) \sim L^{z}$ for both the exponential and integrated autocorrelations, 
where $L^d$ is the lattice size, and $\tau(L)$
is measured at the infinite volume critical point.
For the 2-$d$ Ising model this is known to be 
$\beta_{c} = \log (1 + \sqrt 2) / 2 \approx 0.440687$
\Ising.
For the 3-$d$ model we used the value 0.221654 from ref.
\ref\MCRG3{G.S. Pawley {\it et al.}, Phys. Rev. B {\bf 29} (1984) 4030},
while in the 4-$d$ case we used the value 0.14981 from ref.
\ref\MCRG4{H.W.J. Bl\"ote and R.H. Swendsen, Phys. Rev. B {\bf 22} (1980) 4481}.

Autocorrelations are traditionally measured between each update of the
entire lattice, so for the single cluster Wolff update, where only a fraction
of the lattice sites are updated at each iteration, the autocorrelation time
needs to be scaled by the average cluster size $<|c_{W}|>$ divided 
by the lattice size, so that 
$\tau = \tau' <|c_{W}|> L^{-d}$,
where $\tau'$ is the unscaled autocorrelation time.

The procedure we used was the same as for our previous calculations \us.
We did from 8 to 20 independent runs for each of the different
lattice sizes, with each run being of the order of $10^6$ sweeps.
The results from the multiple runs were averaged to give the overall 
autocorrelation function $\rho_A(t)$, from which the autocorrelation
times were extracted (see ref. \us for details of the fitting procedure).
Each of the individual runs were analyzed in a similar fashion, in order
to provide an error in each of the autocorrelation times. 

For the SW algorithm on the larger lattice sizes in two and three dimensions, 
we used a parallel cluster labeling algorithm which we have developed 
\ref\parallel{C.F. Baillie and P.D. Coddington,
Concurrency: Practice and Experience {\bf 3} (1991) 129}
in order to run on large parallel machines such as the nCUBE hypercube and 
the Symult S2010.
% and the Intel Delta?
For the other lattice sizes, we ran multiple simulations in parallel,
on shared memory machines such as the BBN Butterfly and Encore Multimax, 
or using networks of workstations.

\bigskip
{\bf 3. \quad Results}
\medskip

The results are summarized in Tables 3-5.
For each of the lattice sizes we used, we show
the measured values of $\tau_{int}$ and $\tau_{exp}$ for the 
energy and absolute magnetization, for both the Wolff and SW algorithms,
as well as the total number of sweeps for each algorithm, and the 
measured value of the specific heat.

We found as expected that the energy and absolute magnetization have very
similar values of $\tau_{exp}$, however $|M|$ decorrelates faster than $E$, 
giving smaller values of both $\tau_{int}$ and $z_{int}$ in all cases.

Results for $\tau_{int,E}$ for the different algorithms are plotted in 
Figs. 1,2 and 3 for $d=2$, 3 and 4 respectively.
For $d=3$ we have used a log-log plot, with 
the solid straight lines here representing chi-squared fits to a power law. 
For $d=2$ and $d=4$ we have used a log-linear plot, with 
the solid straight lines representing chi-squared fits to a logarithm. 
Note however that for $d=4$ we plot $\log \tau$ rather than $\tau$ for the 
SW algorithm, so that it will fit on the same graph.
For the cases where alternative behavior is also possible (i.e. a power law
for $d=2$ and $d=4$, and a logarithm for $d=3$), 
we also indicate the alternative fit with dashed lines.
The values of $z$ from fits to all the data are given in Table 6.

In each figure we also include the measured value of the specific heat,
scaled by an arbitrary factor, in order to show that the bound of Li and
Sokal \Li
$$
\tau_{int,E} \ge constant \times C_{H}, z_{int,E} \ge \alpha / \nu
\eqn\Eqbound
$$ 
is indeed satisfied by the Swendsen-Wang algorithm. 
No such bound has been proved for the Wolff algorithm
\ref\Li2{X.-J. Li and A.D. Sokal, Phys. Rev. ???},
although it appears from the figures that not only does the bound hold,
but there may actually be equality in the exponents, so that for the 
Ising model 
$z_{int,E} = \alpha / \nu$ for the Wolff algorithm.

If we compare fits to the specific heat and the Wolff autocorrelation time,
which correspond to the measured values of the exponents $\alpha / \nu$ and
$z_{int,E}$ respectively, then for $d=3$ we get 0.35(1) and 0.33(1).
For $d=2$ and $d=4$ $\alpha = 0$, and $\tau_{int,E}$ for the Wolff 
algorithm can be fitted to a logarithm, implying $z_{int,E}$ is also zero.

It is also apparent from Fig. 1 that the Wolff autocorrelations and the 
specific heat seem to be both increasing logarithmically with $L$, and also
appear to share the same multiplicative constant of the logarithm. The
slope of the fitted lines is 0.64(1) for $C_{H}$ and and 0.61(1) for 
$\tau_{int,E}$.
Note that all the errors quoted are purely statistical (one standard 
deviation), and do not take into account such things as the finite lattice 
size. 
Hence the Wolff algorithm for the 2-$d$ Ising model seems to satisfy the
relation
$$
\tau_{int,E} = constant + C_{H},
\eqn\Eqtau1
$$
while for the 3-$d$ model it appears to satisfy
$$
\tau_{int,E} = constant \times C_{H}.
\eqn\Eqtau2
$$
In four dimensions the relation seems just to be the most general which 
will still allow $z_{int,E} = \alpha / \nu = 0$, which is
$$
\tau_{int,E} = a + b \times C_{H},
\eqn\Eqtau3
$$
where $b$ and $b$ are constants (i.e. independent of $L$). 

Our measurement of $\alpha / \nu$ for the 3-$d$ model gives a value which 
is actually much too large.
This measurement is known to be very difficult
\ref\Landau{D.P. Landau, Phys. Rev. B {\bf 14} (1976) 255}, 
\ref\Barber{M.N. Barber {\it et. al}, Phys. Rev. B {\bf 32} (1985) 1720},
since the divergence is so small that there is a strong
correction to scaling from the non-singular part, which is just the constant
$A$ in the finite-size scaling expression
$$
C_{H} = A + BL^{\alpha / \nu}.
\eqn\EqFSS
$$
Usually when extracting an exponent we would ignore the constant $A$ -- 
including it requires that
we do a three parameter fit to a small number of data points, each of which
has a significant error. Adding this extra parameter generally allows so much
freedom that a wide range of values of the exponent will give a
reasonable fit, so it is very difficult to pin down a precise value for the
exponent from such a fit. 

%
% Try to do this, and compare the improved values of the exponents.
% Is there way to relate the two so that this constant is eliminated?
%

The same problem also occurs in fitting the autocorrelation times to a
power law for the 2-$d$ model. Since the autocorrelations grow so slowly,
it is possible that an additional constant term (as in equation \EqFSS)
could make a big difference, and would certainly improve the fit. 
%
% ******* Try this, and give some results. ***********
% For example, adding ????? and ???? to the values of $\tau_{int,E}$ for the
% Wolff and SW algorithms respectively, reduces the \chi^2/d.o.f. from
% ??? (???) to ??? (???), which matches or betters that of the logarithmic 
% fit. 
% and allows us to fit over a larger range of L???
%
It is therefore very difficult to say with any confidence that $z = 0$, even
with data on very large lattices. The great similarity
between the Wolff autocorrelation times and the specific heat for this model
seems to be the most compelling evidence so far that $z$ is in fact zero
for the 2-$d$ Ising model. 
Since data for the $Q$-state Potts model in two dimensions implies that
$z_{exp}$ and $z_{int,E}$ are the same for the SW and Wolff algorithms
in two dimensions \us, this would imply that $z=0$ for the SW algorithm 
as well.
The equivalence $z_{SW} = z{W}$ 
also implies that the Li-Sokal bound holds for the Wolff algorithm
for 2-$d$ Potts models. However for the case $Q=3$ we find that
$z_{W} \gt \alpha / \nu$ \Li, so the bound is not sharp.

%
% ***************************************************************************
%
% All this stuff needs to be updated.
%
As seen by Burkitt and Heermann, $\tau_{int,|M|}$ seems to
be better fitted with a logarithm, especially for small values of $L$ (although 
these are not really relevant, since we are looking for asymptotic behavior). 
It is therefore quite possible that $z_{int,|M|}$ is zero. 
The other exponents ($z_{int,E}$ and $z_{exp}$ for $E$ and $|M|$) are much 
harder to ascertain, since the autocorrelations
do not give good fits to either a log or a power law unless smaller values of 
$L$ are excluded. 

In Table 4 we give the values of chi-squared per degree of freedom for both 
fits, for the Wolff and SW algorithms. 
The fits include only data with $L \ge 64$ for SW, and $L \ge 50$ for
Wolff, so that we are only fitting to six data points, 
and this may affect the reliability of the chi-squared values. 
The results favor logarithmic behavior, although we would not claim 
that the evidence is conclusive, especially since the differences are so small 
that corrections to scaling could be very important.
% except perhaps for $z_{int}$ of the magnetization for the SW algorithm.
Good data on much larger lattices will probably be necessary to conclusively
differentiate between $z = 0$ and $z \approx 0.25$.
%
% ***************************************************************************
%

\bigskip
{\bf 4. \quad Dynamic Universality}
\medskip

The 3-$d$ Ising model is believed to be in the same static universality 
class as the 3-$d$ 3-state antiferromagnetic Potts model
\ref\???,
\ref\Kotecky{J.-S. Wang, R.H. Swendsen and R. Kotecky, 
Phys. Rev. B {\bf 42} (1990) 2465},
% Phys. Rev. Lett. {\bf 63} (1989) 109},
and the 3-$d$ XY model
\ref\Meyer{M. Hasenbusch and S. Meyer, Phys. Lett. {\bf B241} (1990) 238}.
One would expect that these models would also be in the same {\it dynamic}
universality class, so that each of the cluster algorithms would have
the same dynamic critical exponents for each of the three models.
From Table 7, we see that this appears to be the case for the
SW algorithm, as the exponents are all approximately 0.5, and agree
fairly well within errors (it should be noted that all the errors
quoted are purely statistical, and do not take into account such things
as finite lattice size, so that the actual error may be somewhat larger).
%
%  Ising                0.50(3)  Wolff
%                       0.53(1)  This work
%  XY                   0.46(3)  Hasenbusch and Meyer
%  Potts antiferro      0.48(4)  Wang, Swendsen and Kotecky
%

The case of the Wolff algorithm is particularly interesting, since for
the Ising case, it is uncertain whether the exponent is small, or zero
(logarithmic). 
If we believe that dynamic universality should hold, then we should expect
to get the same dynamical critical exponents for each model, unless there
are non-negligible corrections to scaling, or the exponent is actually 
zero, and what we are seeing is different logarithmic behavior masquerading
as different exponents. Thus different exponents might be evidence for
logarithmic behavior.

We have measured the autocorrelations for the Wolff algorithm for the
3-$d$ antiferromagnetic 3-state Potts model. The results are shown in
Table 9 and Figure 5, where the solid line is a power law fit, and the 
dashed line a fit to a logarithm. We can see that, as in the case of
the Ising model, it is very difficult to distinguish between the two. 
The measured exponents are given in Table ?. 

We also studied the 2-$d$ Ising model on the triangular lattice, to see 
if this model has different dynamical exponents than the square lattice 
model, which may indicate logarithmic behavior. 
The results we obtained were rather surprising -- 
not only were the exponents identical to those of the square lattice model, 
but so were the autocorrelation times,
which are not expected to be universal quantities. In Table ? we show
the autocorrelation times for the triangular lattice Ising model at its
critical temperature of 0.27465307 \Wu. Comparing this to the square 
lattice case (Table 1) shows that the two are identical within errors.

\bigskip
{\bf 5. \quad Conclusions}
\medskip

Klein {\it et al.} \Klein have proposed that 
$$
z_{SW} = z_{G} - 2 \gamma / d_{m} \nu
\eqn\Eqconjecture
$$
where $d_{m}$ is the fractal dimension of clusters, excluding the 
(incipient? largest? percolating? infinite?) cluster, and $z_{G}$ is the
critical exponent for Glauber dynamics. It is difficult to calculate
both $z_{G}$ and $d_{m}$ accurately, however we can assume that $d_{m}$
lies somewhere between the true dimension $d$ and the fractal dimension
of the (incipient? largest? percolating? infinite?) cluster, which is
known to be $d - \beta / \nu$, so that we expect
$$ 
z_{G} - 2 \gamma / d \nu \ge z_{SW} \ge z_{G} - 2 \gamma / (\beta + \gamma)
\eqn\Eqbound
$$
which follows from using the scaling relation $d\nu = 2 \beta + \gamma$.
We can see from the results in Table ?, that the values of $z$ we obtain 
for the SW algorithm satisfy these bounds, unless 
$z_{SW} = 0$ for $d = 2$, or ..................

We can also see that the rigorous bound $z_{SW} \ge \alpha / \nu$ 
obtained by Li and Sokal \Li is satisfied in all cases,
although the bound is far 
from being sharp (with the possible exception of $d = 2$). It is
interesting to note that empirically the Wolff algorithm also appears 
to obey the bound, although this has not been proven. The case $d = 3$ 
could be an exception, since $z_{W} = 0$ is a possibility here. 
Note that if this were a true bound for the Wolff algorithm, it is much
sharper than for SW.

\bigskip
{\bf Acknowledgements}
\medskip

The simulations were done using a number of different computers.
The parallel algorithm was run on a 1024 node Ncube/10, 
and a 192 node Symult S2010, 
while the sequential code was
run on a 16 processor Encore Symmetry, a 20 processor Encore Multimax,
a 96 processor BBN GP1000 Butterfly, 
and a network of SPARCstations, DECstations and IBM RS-6000 workstations.
We would like to thank the North-East Parallel Architecture Center at
Syracuse University, 
Michigan State University, 
the Caltech Concurrent Supercomputer Facility,
and the Sandia National Laboratories for the use of these machines, 
which were accessed under the auspices of the Center for Research in 
Parallel Computing.
This work was sponsored in part by Department of Energy grants
DE-FG03-85ER25009 and DE-AC03-81ER40050, 
and by a grant from the IBM corporation.
CFB is supported by Air Force Office of Scientific Research Grant AFOSR-89-0422
and by the Department of Energy under contract DE-AC02-86ER40253.

\listrefs\vfill
\end
