\input /home/ace/paulc/texlib/mac.tex
\input /home/ace/paulc/texlib/moremac.tex
\input /home/ace/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-TR91193 \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
% For double space letter
% \smallskip
\centerline{\it Physics Department, Syracuse University, Syracuse NY 13244, U.S.A.}
\bigskip
\bigskip
\centerline{\big Clive F. Baillie} 
% \medskip
% For double space letter
% \smallskip
\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 dynamic 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
% For double space letter
% \vskip 2cm
\noindent
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 (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 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)}).
%
There is little theoretical understanding of the dynamics of these 
algorithms. In particular, little is known as to why they seem to 
eliminate critical slowing down completely in some cases,
and not others. There is no known theory which can predict the value of 
the dynamic critical exponent $z$ for any spin model, although a 
rigorous bound on $z$ for the SW algorithm for Potts models 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)}.

The measurement of dynamic critical exponents is notoriously 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. 

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 dynamic
critical exponents, we have done numerical 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 dynamic exponents.
We measured the time correlation function $\rho(t)$ for the energy, 
and extracted the integrated autocorrelation time
\Sokalrev
$ \quad \tau = {1 \over 2} + \sum_{t=1}^{\infty} \rho(t)$.
The dynamic critical exponent $z$ is given by $\tau \sim L^z$, 
where $\tau$ for the different lattice sizes is measured 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)}.
A detailed account of the methods we used to do the measurements, fits and
error estimates, is given in ref. \us.

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
$$
is what we present for the Wolff autocorrelations. 
Since this scaling ratio is an estimator for the susceptibility \Wolff,
the dynamic critical exponent $z'$ for the unscaled autocorrelations is 
given by ${z' = z + (d - \gamma / \nu)}$, where 
$\nu$ is the critical exponent for the correlation length,
and $\gamma$ is the critical exponent for the susceptibility, which 
diverges as $L^{\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.
For the other lattice sizes, we ran multiple simulations in parallel
using smaller shared memory machines and 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 very difficult 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. Here $\alpha$ is the critical
exponent for the specific heat ($C_H \sim L^{\alpha/\nu}$).
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.

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$, 
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. 
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).
In two dimensions the best fit is obtained with ${a \approx -0.474}$ and 
${b \approx 0.957}$ (the data does not exclude the possibility that 
${b = 1}$, which would imply that $\tau_{int,E}^{W}$ is just a constant
plus $C_{H}$). 
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 exponent $\beta / \nu$ characterizing the divergence of the
magnetization has values which are similar
to our measured values for the dynamic 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 
constant, 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), and fit very poorly to a power law.
The data therefore support 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
$$

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 dynamic exponents, as expected from 
relations \EqWtau\ and \EqSWtau.
Thus although we may not be able to measure the asymptotic behavior
of the autocorrelation times, finding simple relations between
the autocorrelations and static quantities whose asymptotic behavior
is known enables us to infer the true values of the dynamic critical
exponents.

This is especially useful for the 2-$d$ model, for which 
the autocorrelations grow so slowly that any corrections to scaling 
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\ 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
of the SW and Wolff cluster algorithms for the Ising model in 2, 3 and
4 dimensions. We have found what appear to be surprisingly simple empirical 
relations between the autocorrelation times of these algorithms and 
simple static quantities (the magnetization and specific heat). 
These relations could perhaps stem from the fact that the dynamics of 
cluster algorithms are closely linked to 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 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.
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}.

\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 constant. 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
