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

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

% *** for double-spaced letter ***
% \baselineskip=24pt
%

% *** for preprint ***
\line{SCCS-141 \hfill}
% \line{CRPC-TR91??? \hfill}
\line{COLO-HEP-264 \hfill}
\bigskip
\bigskip
%
\bigskip
\centerline{\bfbigger Empirical relations between static and dynamic} 
\smallskip
\centerline{\bfbigger exponents for Ising model cluster algorithms} 
\bigskip
\bigskip
\bigskip
\centerline{\big Paul D. Coddington}
\medskip
\centerline{\it Physics Department, Syracuse University, Syracuse NY 13244, U.S.A.}
\bigskip
\bigskip
\centerline{\big Clive F. Baillie} 
\medskip
\centerline{\it Physics Department, University of Colorado, Boulder CO 80309, U.S.A.}
\bigskip
\bigskip
\bigskip
% *** for preprint ***
\centerline{September 1991}
\bigskip
\bigskip
\bigskip
%
\bigskip
\centerline{\bf Abstract}
\medskip
\noindent
We have measured the autocorrelations 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 the autocorrelations are
linearly related to the specific heat, in which case the dynamical critical 
exponent ${z_{int,E}^{W} = \alpha / \nu}$.
For the Swendsen-Wang algorithm, scaling the autocorrelations by the
average maximum cluster size gives either a constant or a logarithm,
which implies that ${z_{int,E}^{SW} = \beta / \nu}$ for the Ising model.

% *** for letter ***
% \vskip 5cm
% PACS numbers 05.50.+q, 11.15.Ha, 64.60.Ht.
%

\vfill
\eject


\bigskip
{\bf 1. 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}, 86 (1987)}
and Wolff
\ref\Wolff{U. Wolff, Phys. Rev. Lett. {\bf 62}, 361 (1989)}
can dramatically reduce critical slowing down in computer simulations of spin 
models, and thus greatly increase the computational efficiency of the
simulations. 
Since they were introduced, a large amount of effort has gone into
studying these algorithms, in an attempt to understand how they work, 
and how they can be generalized to more complicated models
(for reviews of cluster algorithms, see refs.
\ref\Sokalrev{A.D. Sokal, in: {\it 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}, 55 (1991)}
\ref\reviews{U. Wolff, Proc. of the 1989 Symposium on Lattice Field
Theory, Capri, (Sept. 1989), Nucl. Phys. B (Proc. Suppl.) {\bf 17}, 93 (1990);
J.-S. Wang and R.H. Swendsen, Physica A {\bf 167}, 565 (1990);
C.F. Baillie, Int. J. of Mod. Phys. C {\bf 1}, 91 (1990)}).

So far there has been very little theoretical understanding of the dynamics
of these algorithms. In particular, little is known as to why cluster
algorithms seem to eliminate critical slowing down completely in some cases,
and not in others. There is no known quantitative theory which can predict
the value of the dynamical critical exponent $z$ for any spin model, although
a rigorous bound on $z$ for the SW algorithm has been derived
\ref\Li{X.-J. Li and A.D. Sokal, Phys. Rev. Lett. {\bf 63}, 827 (1989)}.
Another problem which is not well understood is why the SW and Wolff 
algorithms give similar values of $z$ for the 2-$d$ Potts model
\ref\us{C.F. Baillie and P.D. Coddington, Phys. Rev. B {\bf 43}, 10617 (1991)},
but have very different behavior for other models, such as the Ising
model in more than two dimensions 
\ref\WIsing{U. Wolff, Phys. Lett. {\bf B228}, 379 (1989)}
\ref\Tamayo{P. Tamayo, R.C. Brower and W. Klein,
J. Stat. Phys. {\bf 58}, 1083 (1990)}.

A lot of work has been done on measuring the dynamical critical exponents
of cluster algorithms for various spin models. These results are very useful
``experimental'' inputs for any theoretical understanding of these algorithms,
since $z$ is a quantity which one would hope to be able to derive, at least 
approximately, from any comprehensive theory of cluster algorithm dynamics.
Unfortunately the measurement of dynamical critical exponents is very 
difficult, and both very good statistics and very large lattices are required
in order to obtain accurate results. 
This is certainly the case for the Ising model, where a number of different
measurements have given conflicting results. 

For the two dimensional Ising model, initial results suggested
${z \approx 1/3}$ for both the SW and Wolff algorithms \SW \Tamayo. 
Further work \WIsing\ gave ${z \approx 1/4}$, and it was later shown that 
the data were consistent with a logarithmic divergence, suggesting that 
${z = 0}$
\ref\Burkitt{D.W. Heermann and A.N. Burkitt, Physica A {\bf 162}, 210 (1990)}.
Recent results show that it is very difficult to distinguish between
a logarithm and a small power \us.
Measurements on the three dimensional model have proven to be just as 
difficult, with values of $z$ for the SW algorithm ranging from 
0.339(4) to 0.75(1) \SW \WIsing 
\ref\Wang{J.-S. Wang, Physica A {\bf 164}, 240 (1990)}.
For the Wolff algorithm, Tamayo {\it et al.} \Tamayo\ obtained 0.44(10), 
while Wolff found a value of 0.28(2) for the energy autocorrelations \WIsing.
We have examined Wolff's data and found that it also fits well to a logarithm,
so that ${z = 0}$ is also a possibility. 
Again it is very hard to tell the difference between a logarithm
and a small power.
In four dimensions only one result is 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
expected to give the same exponents as the Ising model in four or more
dimensions
\ref\Ma{S.-K. Ma, {\it Modern Theory of Critical Phenomena},
(Benjamin, New York, 1976)}.
The mean-field data are consistent with $z$ being 0 for the Wolff algorithm
\Tamayo\ and 1 for SW 
\ref\MF{T.S. Ray, P. Tamayo and W. Klein, Phys. Rev. A {\bf 39}, 5949 (1989)},
with the latter result being supported by theoretical arguments.

\bigskip
{\bf 2. Simulations}
\medskip

Due to the discrepancies between the various measurements of the dynamical
critical exponents, 
we have been performing simulations of the Ising model in
2, 3 and 4 dimensions using the SW and Wolff algorithms,
with the aim of obtaining good statistics on fairly large lattices, 
in order to get reliable values for the dynamical critical exponents.
The definitions of the various autocorrelations and dynamical exponents, 
and an account of the methods we used to do the measurements, fits and
error estimates, are given in ref. \us.
We use finite size scaling to determine $z$, by measuring the autocorrelations
for different lattice sizes at the infinite volume critical point.
We have used the Potts formulation of the Ising model, 
for which the critical point in two dimensions is known to be 
$\beta_{c} = \log (1 + \sqrt 2) \approx 0.8813736$
\ref\Wu{F.Y. Wu, Rev. Mod. Phys. {\bf 54}, 235 (1982)}.
For the 3-$d$ model we used the value 0.443308 
\ref\MCRGthree{G.S. Pawley {\it et al.}, Phys. Rev. B {\bf 29}, 4030 (1984)},
while in the 4-$d$ case we have used 0.29962 
\ref\MCRGfour{H.W.J. Bl\"ote and R.H. Swendsen, 
Phys. Rev. B {\bf 22}, 4481 (1980)}.

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 measured
autocorrelation time $\tau'$
needs to be scaled by the ratio of the average Wolff cluster size $<|c_{W}|>$ 
and the number of lattice sites $L^d$. The scaled autocorrelation time
$$
\tau = \tau' <|c_{W}|> / L^d
\eqn\Eqscale
$$
gives a more reasonable comparison between the Wolff and SW algorithms in
terms of the amount of work done (or CPU time used), and it is these scaled 
autocorrelations that we present here for the Wolff algorithm.
Since this scaling ratio is an estimator for the susceptibility
\Wolff,
the dynamical critical exponent $z'$ for the unscaled autocorrelations is 
given by ${z' = z + (d - \gamma / \nu)}$.

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}, 129 (1991)}
in order to run on large parallel machines such as the nCUBE hypercube and 
the Symult S2010.
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. Results}
\medskip

Results for $\tau_{int,E}$, the integrated autocorrelation time for the energy,
are shown in Figures 1(a), (b) and (c) for $d=2$, 3 and 4 respectively.
For $d=3$ we have used a log-log plot, with 
the straight lines representing $\chi^2$ fits to a power law, while
for $d=2$ and $d=4$ we have used a log-linear plot, with 
the straight lines representing $\chi^2$ fits to a logarithm. 
Note however that for $d=4$ we plot $\log \tau_{int,E}$ rather than 
$\tau_{int,E}$ for the SW algorithm, since the SW autocorrelations 
increase as a power of $L$.
The measured values of the exponents from the fits to the data are shown
in Table 1.
For the Wolff algorithm in all dimensions, and the SW algorithm in two 
dimensions, it is virtually impossible to distinguish between a small
exponent and a logarithmic increase in the autocorrelations 
(which would imply that $z = 0$). 

% *** for preprint ***
\bigskip
\bigskip
\input tab1.tex
\bigskip
\bigskip
% 
In Fig. 1 we also include the measured value of the specific heat
$C_H$, scaled by an appropriate factor, in order to show that the bound of 
Li and Sokal \Li
$$
\tau_{int,E} \ge constant \times C_{H}, \qquad z_{int,E} \ge \alpha / \nu
\eqn\Eqbound
$$ 
is indeed satisfied by the SW algorithm. 
No such bound has been proven for the Wolff algorithm
\ref\LiSok{X.-J. Li and A.D. Sokal, private communications},
although it appears from the figures that not only does the bound hold,
but that there may actually be equality in the exponents, so that for the 
Ising model in all dimensions
${z_{int,E} = \alpha / \nu}$ for the Wolff algorithm.

If we compare the results of fits to $C_H$ and $\tau_{int,E}^{W}$
(the autocorrelations in the energy for the Wolff algorithm),
which correspond to the measured values of 
$\alpha / \nu$ and $z_{int,E}^{W}$ respectively, then for 
$d=3$ we find 0.32(1) and 0.33(1).
In two and four dimensions $\alpha = 0$, 
% \Ma, \???
and $z_{int,E}^{W}$ is also consistent with zero.
Hence the Wolff algorithm for the Ising model seems to satisfy the
surprisingly simple relations
$$
\tau_{int,E}^{W} = a + b \times C_{H}, \qquad z_{int,E}^{W} = \alpha / \nu,
\eqn\EqWtau
$$
where $a$ and $b$ are constants (i.e. independent of $L$). 
In Figure 2 we plot the difference ${\tau_{int,E}^{W} - (a + b \times C_{H})}$
for the various dimensions, with $a$ and $b$ chosen to minimize
$\chi^2$ over a certain range of lattice sizes (smaller values of $L$ are
excluded from the fit). We can see that in all cases, values of $a$ and $b$ 
can be found such that the difference is zero within errors.
Note that all the errors shown here are purely statistical (one standard
deviation).

It appears from Fig. 1(a) that in 2-$d$ the Wolff autocorrelations and the 
specific heat share the same multiplicative constant of the 
logarithm, since the lines of best fit have approximately the same slope. 
This corresponds to $b=1$, in which case $\tau_{int,E}^{W}$ is just a constant
plus $C_{H}$. From our data this is certainly a possibility, however
the best fit is obtained with ${a \approx -0.474}$ and ${b \approx 0.957}$, 
and this is what we show in Fig. 2(a).
For the 3-$d$ model the additive constant $a$ is consistent with zero, so that
the autocorrelation time may be just a multiple of the specific heat, with
${b \approx 0.148}$.
In four dimensions we find ${a \approx 0.167}$ and ${b \approx 0.050}$.

The surprising simplicity of the result \EqWtau\ led us to look for a similar
relation for the SW algorithm. 
The power of cluster update algorithms comes from the fact that they flip 
large clusters of spins at a time.
The relative average size of the largest SW cluster, 
${m = <|c_{SW}^{max}|> / L^d}$, is an estimator of the magnetization
\ref\Binder{M. D'Onorio De Meo, D.W. Heermann and K. Binder,
J. Stat. Phys. {\bf 60}, 585 (1990)},
and the magnetization exponent $\beta / \nu$ has values which are similar
to our measured values for the dynamical exponents of the SW algorithm.
To investigate this further, we have scaled the SW autocorrelations by $m$,
in a similar manner to the scaling of the Wolff autocorrelations in equation
\Eqscale. If this gives a constant or a logarithm, then $\tau_{int,E}^{SW}$
diverges like the magnetization, and so we have 
${z_{int,E}^{SW} = \beta / \nu}$.

The SW autocorrelations scaled by $m$ (and by an additional arbitrary factor,
so that these points are not entangled with others in the plots) are also
shown in Fig. 1. 
For ${d = 4}$ the results are very close to a constant, while for ${d = 3}$ they
seem to approach a constant as $L$ increases. 
In two dimensions the scaled autocorrelations are not constant, but they
fit very much better to a logarithm than does the unscaled data, 
as can be seen in Fig. 1(a).
The data therefore supports the assertion that
$$
m \hskip 0.5em\relax \tau_{int,E}^{SW} = a + b \times \log L, \qquad z_{int,E}^{SW} = \beta / \nu.
\eqn\EqSWtau
$$
We are still faced with the problem that it is difficult to 
tell if the scaled autocorrelations increase as a logarithm or a small power,
so that unless ${b = 0}$ we can only say with any certainty that 
${z_{int,E}^{SW} \approx \beta / \nu}$. 
However the scaled 2-$d$ data fit very much better to a logarithm
than to a power law.

Our measurements of $z_{int,E}^{W}$ in 3-$d$ and
$z_{int,E}^{SW}$ in 4-$d$ give results which are very different from
the accepted values of $\alpha / \nu \hskip 0.5em\relax (\approx 0.10 - 0.20)$ 
\ref\Landau{D.P. Landau, Phys. Rev. B {\bf 14}, 255 (1976)}
\ref\Barber{M.N. Barber {\it et. al}, Phys. Rev. B {\bf 32}, 1720 (1985)}
and $\beta / \nu \hskip 0.5em\relax (= 1)$
\MCRGfour
\ref\Mouritsen{O.M. Mouritsen and G.J. Knak Jensen, Phys. Rev B {\bf 19},
3663 (1979)},
since corrections to scaling are known to be important for these
quantities. If we do a simple power law fit to the specific heat in
3-$d$ and the magnetization in 4-$d$, we get results which are also very
different from the actual exponents, but which are very close to the 
measured values of the corresponding dynamical exponents, as expected from 
relations \EqWtau\ and \EqSWtau.
Thus although we may not be able to measure the asymptotic behavior
of the autocorrelation times, if we can find simple relations between
the autocorrelations and static quantities whose asymptotic behavior
is known, then we can infer the true value of the dynamical critical
exponents.

This is especially useful for the 2-$d$ model, for which 
the autocorrelations grow so slowly that any corrections to scaling 
% (such as a simple additive constant) 
could have a big effect.
It is therefore very difficult to say with any confidence that $z = 0$, even
with data on very large lattices. The apparent relation \EqWtau\
between the Wolff autocorrelation times and the specific heat seems to be 
the most compelling evidence so far that $z_{int,E}^{W}$ is in fact zero
for the 2-$d$ Ising model, while the relation \EqSWtau\ would imply that
$z_{int,E}^{SW}$ is actually 1/8, which is not apparent from the
usual fits to either a logarithm or a power law.

\bigskip
{\bf 4. Conclusions}
\medskip

We have measured the autocorrelations and dynamic critical exponents
for the SW and Wolff cluster algorithms for the Ising model in 2, 3 and
4 dimensions. We have found what appear to be very simple empirical 
relations between the autocorrelation times of these algorithms and 
simple static quantities (the magnetization and specific heat). 

The simplicity of the relations \EqWtau\ and \EqSWtau\ between static and 
dynamic quantities is very surprising. 
It could perhaps be due to the fact that the dynamics of these algorithms
are closely linked with the physical properties of the system, 
since the Swendsen-Wang clusters are just the Coniglio-Klein-Fisher droplets
\ref\Fisher{M.E. Fisher, Physics {\bf 3}, 255 (1967);
A. Coniglio and W. Klein, J. Phys. A {\bf 13}, 2775 (1980)},
or ``physical clusters'' \Binder,
from which the critical behavior of the system may be described.
The simplicity of these relations indicates that a quantitative theoretical
understanding of the dynamic critical behavior of cluster algorithms in 
terms of fundamental static quantities should be attainable, at least for 
the Ising model.
It is of course possible that these relations are not exact, but merely
good approximations.
We are currently collecting more data in order to check whether these
results hold up with larger lattices and better statistics, and we will 
present more detailed results in a future publication
\ref\next{P.D. Coddington and C.F. Baillie, in preparation}.

The relations \EqWtau\ and \EqSWtau\ are certainly not general results, 
since for the 2-$d$ $q=3$ Potts model we find that
${z^{W} > \alpha / \nu}$ and ${z^{SW} > \beta / \nu}$ 
\Li \us. 
Also, it is quite surprising that these empirical relations imply that 
$z^{SW}$ is not equal to $z^{W}$ for the 2-$d$ Ising model, whereas
the two appear to be equal for the 2-$d$ $q=3$ Potts model.
However it is interesting to note that for this model
${(\alpha + \beta) / \nu = 8/15 \approx 0.5333}$, which is very close to
the measured values for $z_{int,E}$, of 0.55(1) for SW and 0.57(1) for
Wolff \us.

It would be interesting to investigate whether simple relations between the 
autocorrelations of cluster algorithms and static quantities such as the
specific heat and the magnetization can also be found for other spin models, 
such as Potts and $O(N)$ models. 

\bigskip
{\bf Acknowledgements}
\medskip

The simulations were done using an nCUBE/10, a Symult S2010, 
an Encore Multimax, a 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, 
the Caltech Concurrent Supercomputer Facility,
Michigan State University, 
and the Sandia National Laboratories for the use of these machines,
as well as the Center for Research in Parallel Computing for allowing
us access to most of these computers.
PDC 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.

We would also like to thank Alan Sokal for useful discussions.

\listrefs
\vfill\eject

% *** for letter ***
% \input tab1.tex
% \vfill\eject
%

% Figure captions
\centerline{Figure captions}
\bigskip
\bigskip

\noindent
Fig. 1. \quad Autocorrelations $\tau_{int,E}$ for the Wolff and SW algorithms 
plotted against lattice size $L$ for the Ising model in
(a) 2-$d$, (b) 3-$d$ and (c) 4-$d$.  
Also shown is the specific heat $C_H$, and the SW autocorrelations
scaled by the average maximum cluster size $m$. The latter two quantities
are also scaled by an arbitrary factor. The plots are log-linear for (a) and
(c), and log-log for (b). All error bars are shown, but are usually smaller
than the points. The lines are fits to a power law, logarithm, or constant.

\bigskip
\bigskip

\noindent
Fig. 2. \quad The difference $\tau_{int,E}^W - (a + b \times C_H)$ between
the Wolff autocorrelations and a simple linear function of the specific
heat, for the Ising model in (a) 2-$d$, (b) 3-$d$ and (c) 4-$d$.
The values of $a$ and $b$ are chosen so as to minimize the $\chi^2$, except
in three dimensions, where we have taken $b = 0$. The errors shown are almost
all less than 1\% of $\tau_{int,E}^W$.

\end
