\input   mac.tex
\input   moremac.tex
\input peqns.tex

\line{\bf 1. Introduction\hfill}

Spin models were invented as simple statistical mechanical models of
ferromagnetism. In most cases they exhibit the cooperative behavior found in
phase transitions, which arises from the development of long range order in the system.
This occurs at a certain critical temperature called the Curie point.
The simplest example of such a model, in which the magnetic moments are assumed to be
classical one-dimensional spins capable of only two orientations,
is the so-called Ising model
\ref\Ising{E. Ising, {\it Z. Physik} {\bf 31}, 253 (1925)}
which was actually proposed by Ising's research director Lenz
\ref\Lenz{W. Lenz, {\it Physik Z.} {\bf 21}, 613 (1920)}.
Ising carried out an exact calculation of the partition function for this model
in {\it one dimension} and showed that there was no phase transition to a
ferromagnetic state at any temperature.
Heisenberg
\ref\Heisenberg{W. Heisenberg, {\it Z. Physik} {\bf 49}, 619 (1928)}
used this supposed failure of the Ising model to explain ferromagnetism as one
justification for developing his own theory based on a more complicated model, in
which the magnetic moments are related to quantum-mechanical three-component spin
operators.  Finally, Onsager
\ref\Onsager{L. Onsager, {\it Phys. Rev.} {\bf 65}, 117 (1944)}
solved the {\it two-dimensional} Ising model and demonstrated that there was indeed
a ferromagnetic phase transition.
We now believe that both the Ising model in greater than one dimension and the
Heisenberg model in greater than two dimensions have second order phase transitions.
This belief is based mainly on numerical computer simulations since,
unfortunately, there are as yet no analytical solutions for the Ising model in more
than two dimensions or the Heisenberg model in more than one dimension.

In this review, we are going to focus on computer simulations of
spin models. In these, the spins are usually set up on the sites of a $d$-dimensional
hypercubic lattice of length $L$. The $L^d$ spins form some configuration.
The goal of computer simulations is to generate spin configurations typical of statistical
equilibrium and measure physical observables on this ensemble of configurations.
We shall restrict ourselves to looking at the generation part, since recently there
has been a great deal of activity in producing efficient new algorithms for this.
The main obstacle to efficiency is the phenomenon known as critical slowing down.

Generation of configurations is traditionally performed by Monte Carlo methods.
The first application of these to problems in statistical mechanics was by
Metropolis {\it et al} \ref\Metro{
N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller
and E. Teller, {\it J. Chem. Phys.} {\bf 21}, 1087 (1953)}.
All Monte Carlo work has the same general structure: given some probability
distribution $p$, we wish to generate many random samples $\phi$ from $p$.
Monte Carlo does this by constructing a Markov chain having $p$ as its
equilibrium distribution. In other words, we invent a transition probability
$p(\phi \rightarrow \phi')$ that leaves $p$ invariant
$$\int dp(\phi) p(\phi \rightarrow \phi') = dp(\phi') \eqno\eq$$
and also satisfies an ergodicity condition
(that is, it is possible to get from any state to any other).
A sufficient condition for this Markov process to converge to $p$ is that the
transition probabilities obey detailed balance
$$dp(\phi) p(\phi \rightarrow \phi') = dp(\phi') p(\phi' \rightarrow \phi) \eqno\eq$$
(that is, the probability of going from one state to another is the same in
both directions).
In physical terms, we are inventing a stochastic time evolution for the given system
(which is actually very similar to the evolution produced by Langevin dynamics
\ref\MetLanT{P. Meakin, H. Metiu, R. G. Petschek and D. J. Scalapino, {\it J. Chem. Phys.}
{\bf 79}, 1948 (1983); C. F. Baillie and D. A. Johnston, \PRD{39}, 1246 (1989)} ).
Note, however, that the Monte Carlo time evolution need not correspond to any real
physical dynamics; rather, the dynamics is simply a numerical algorithm and therefore
is chosen on the basis of its computational efficiency.
For a statistical mechanical system such as a spin model the probability distribution
we require is the Boltzmann distribution, hence $p = e^{- \beta S(\phi)}$, where
$S(\phi)$ is the action, or energy, of the system in configuration $\phi$,
and $\beta = 1/k_B T$ is the inverse temperature, or coupling (in numerical computations
we shall set the Boltzmann constant $k_B = 1$).
Thus, on the computer, we set up an arbitrary initial configuration $\phi_0$ and use
Monte Carlo to generate the Markov chain $\phi_1, \phi_2, ...$. After some
equilibration time of $E$ steps, we calculate the expectation value of the operator $O$
representing a physical observable from the next $T$ configurations as
$$<O> = {1 \over T} \sum_{i=E+1}^{E+T} O(\phi_i). \eqn\EqM$$
Markov chain theory
\ref\Markov{J. G. Kemeny and J. L. Snell, {\it Finite Markov Chains} (Springer, New York,
1976); M. Iosifescu, {\it Finite Markov Processes and Their Applications} (Wiley,
Chichester, 1980)} guarantees that such time averages converge to $p$-averages so that
\EqM\ is equivalent to the functional integral we want
$$<O> = {1 \over Z} \int D\phi O(\phi) e^{- \beta S(\phi)}, \eqno\eq$$
where the partition function
$$Z = \int D\phi e^{- \beta S(\phi)}. \eqno\eq$$
The trouble in practice is that the samples $\phi_i, \phi_{i+1}, ...$ from the Markov
chain are not statistically independent but are rather correlated.
The statistical error in Monte Carlo behaves as $1/\sqrt{N}$, where $N$ is the number of
{\it effectively independent} samples. Thus the statistical error is $1/\sqrt{T/2\tau}$, where
$\tau$ is the autocorrelation time of the Markov chain.
For a spin model with a phase transition, as the inverse temperature $\beta$ approaches
the critical inverse temperature $\beta_c$, $\tau$ diverges to infinity -- so that the
computational efficiency goes to zero! This behavior is called critical slowing down (CSD).
Until very recently this problem has plagued Monte Carlo simulations of statistical
mechanical systems, in particular spin models, at or near their phase transitions.
Now, however, there are several new algorithms which manage to avoid CSD
-- partially or even completely -- thus facilitating better computer simulations.
We shall describe these algorithms in detail in Secs. 4 and 5
(explaining how they get around CSD, whereas
the older Monte Carlo methods do not); after setting the scene with a brief historical
review of the numerical computer work performed on spin models in Sec. 2,
including a summary of the measurements of CSD in Sec. 3.

\bigskip
\line{\bf 2. History\hfill}

The partition function for a general spin model, in $d$ dimensions, may be written
$$Z = \int D\sigma e^{- \beta S(\sigma)}, \eqno\eq$$
where the action
$$S = -J \sum_{<ij>} s(\sigma_i, \sigma_j), \eqno\eq$$
with $s(\sigma_i, \sigma_j)$
being the energy of a nearest neighbor pair,
denoted $<ij>$, of spins on the lattice.
If the spins are pointing in the same direction then this energy is $-J$ so
the system is ferromagnetic for $J > 0$ (we set $J = 1$ on the computer).
The sum is restricted to nearest neighbor
interactions for computational simplicity -- this restriction does not affect the
critical behavior of the model, that is, the behavior near a phase transition (PT).
There are three main features of PTs in which we are interested:
\medskip
\item{i.} where it is, i.e., at what critical temperature $T_c$
\smallskip
\item{ii.} what order it is, i.e., first, second or Kosterlitz-Thouless
\smallskip
\item{iii.} what the critical exponents are, i.e., what exponents do physical
observables diverge with at the critical temperature
\medskip
\noindent
(The reader is assumed to be familiar with this jargon.)
The $\sigma_i$'s may be discrete or continuously valued variables depending on the model.

For q-state Potts models
\ref\Potts{R. B. Potts, {\it Proc. Camb. Phil. Soc.} {\bf 48}, 106 (1952)},
$s(\sigma_i, \sigma_j) = \delta_{\sigma_i, \sigma_j}$ with $\sigma_i$ taking one
of $q$ discrete values $0,1,...,q-1$.
For $q=2$ the Potts model is equivalent to the Ising model whose action is usually written
$$S = -J \sum_{<ij>} \sigma_i \sigma_j, \eqn\EqIsing$$
with the spins taking on the values $+1$ (up) or $-1$ (down).
Historically, a four-component version of the model was first studied by Ashkin and Teller
\ref\AaT{J. Ashkin and E. Teller, {\it Phys. Rev.} {\bf 64}, 178 (1943)}. However, the
general q-component model bears its current name after it was proposed by Domb
\ref\Domb{C. Domb, {\it J. Phys. A} {\bf 7}, 1335 (1974)}
to his research student Potts as a thesis topic.
The phase transitions exhibited by the q-state Potts model in various dimensions are
shown in Table 1.
\medskip
\vbox{
\bigskip
\input sptab1
\bigskip
}

Alternatively, if the spins are unit vectors in $N$-dimensions
($\sigma_i \in S_{N-1}$)
then we obtain the O(N) models with $s = \sigma_i . \sigma_j$.
$N = 1$ gives the Ising model again. $N = 2$ is the XY model, also known as the planar
Heisenberg or planar rotator model (first worked on by Vaks and Larkin
\ref\VaL{V. G. Vaks and A. I. Larkin, {\it Zh. eksp. teor. Fiz.} {\bf 49}, 975 (1965)
[translation: {\it Soviet Phys. JETP} {\bf 22}, 678 (1966)]}),
in which the dot product reduces to the cosine of the
angle between the spins which are simply represented as angles: $s =
cos(\theta_i - \theta_j)$, with $\sigma_i = e^{i\theta_i}$.
$N = 3$ is the classical Heisenberg model \Heisenberg.
Finally, as $N \rightarrow \infty$, we get the spherical model
\ref\Sphere{T. H. Berlin and M. Kac, {\it Phys. Rev.} {\bf 86}, 821 (1952)}.
The phase transitions exhibited by the O(N) Heisenberg model in various dimensions are
shown in Table 2.
\medskip
\vbox{
\bigskip
\input sptab2
\bigskip
}

(There is another class of models, called the planar Potts, vector Potts or clock model
-- denoted Z(q) --
in which the spins are confined in a plane with each spin pointing to one of the $q$
equally spaced directions specified by the angles $\theta_n = 2 \pi n / q, n=0,1,...,q-1$,
and with $s = cos(\theta_{n_i} - \theta_{n_j})$. We are not going to discuss
these (see
\ref\Wu{F. Y. Wu, {\it Rev. Mod. Phys.} {\bf 54}, 235 (1982)}
for a review),
other than noting that Z(2) is of course the Ising model, and Z($\infty$) is the XY model.)

\bigskip
\line{\bf 2.1. q-state Potts models\hfill}

Most important information about the Potts model can be found in the excellent review
by Wu \Wu.
He points out, for example, that the critical point is now rigorously known for
two-dimensional square, triangular and honeycomb lattices for all $q \geq 4$
\ref\Hinter{A. Hintermann, H. Kunz and F. Y. Wu, {\it J. Stat. Phys.} {\bf 19}, 623 (1978)},
as well as for $q = 2$ \Onsager\ of course.
The critical condition for $q \geq 4$ is expected to hold for $q = 3$ (since it agrees
for the exact $q = 2$ case) but this has not been proven. 
Baxter \ref\Baxter{R. J. Baxter, {\it J. Phys. C} {\bf 6}, L445 (1973)}
has shown that this transition is first order for $q > 4$ and second order for $q \leq 4$.
This has been verified by the numerical calculations of Rebbi and Swendsen
\ref\RaS{C. Rebbi and R. H. Swendsen, \PRB{21}, 4094 (1980)} for $q = 3,4,6,7$, and Binder
\ref\Binder{K. Binder, {\it J. Stat. Phys.} {\bf 24}, 69 (1981)} for $q = 3,4,5,6$.
Later, Swendsen, Andelman and Berker
\ref\SAN{R. H. Swendsen, D. Andelman and A. N. Berker, \PRB{24}, 6732 (1981)} looked at
$q = 4$ in more detail, and Swendsen and Berker
\ref\SN{R. H. Swendsen and A. N. Berker, \PRB{28}, 3897 (1983)} carefully re-examined
$q = 3$.
Finally conjectures
\ref\conjs{M. P. M. den Nijs, {\it J. Phys. A} {\bf 12}, 1857 (1979);
B. Nienhuis, A. N. Berker, E. K. Riedel and M. Schick, \PRL{43}, 737 (1979);
B. Nienhuis, E. K. Riedel and M. Schick, {\it J. Phys. A} {\bf 13}, L189 (1980);
R. B. Pearson, \PRB{22}, 2579 (1980)}
for the exact critical exponents of all two-dimensional Potts models have been proven
\ref\proofs{J. L. Black and V. J. Emery, \PRB{23}, 429 (1980);
B. Nienhuis, {\it J. Phys. A} {\bf 15}, 199 (1982);
M. P. M. den Nijs, \PRB{27}, 1674 (1983)};
a good explanation of all this is
\ref\gooexp{A. N. Berker and D. Andelman, {\it J. Appl. Phys.} {\bf 53}, 7923 (1982)}.

In three or more dimensions no exact results are known and so the critical temperature must
be located by numerical means.
The 3-state Potts model is the most interesting since in 3d since the order
of the phase transition is not clear (see Table 1).
One of the first three-dimensional simulations of the 3-state Potts model was performed by
Herrmann \ref\Herrmann{H. J. Herrmann, {\it Z. Physik} {\bf B 35}, 171 (1979)}
in 1979 (on a CDC CYBER 76 computer); using lattices up to $20^3$ he found a first order
PT at $T_c = 1.818 \pm 0.003$.
Around the same time, Bl\"ote and Swendsen
\ref\BaS{H. W. J. Bl\"ote and R. H. Swendsen, \PRL{43}, 799 (1979)} simulated the same
model in three and four dimensions, obtaining a first order PT at
$T_c = 2.581 \pm 0.007$ in 4d.
Moreover, Knak Jensen and Mouritsen
\ref\KJaM{S. J. Knak Jensen and O. G. Mouritsen, \PRL{43}, 1736 (1979)}
also concurrently found a first order PT for the 3-state model, with much higher statistics
using a larger lattice ($30^3$).
More recently, Gavai, Karsch and Petersson
\ref\GKaP{R. V. Gavai, F. Karsch and B. Petersson, ``A Study of the Correlation Length
near a First-Order Phase Transition: The 3-d Three-State Potts Model'', CERN-TH-5221/88
Preprint (November 1988)} concur.
On the other hand, Fucito and Vulpiani
\ref\FaV{F. Fucito and A. Vulpiani, {\it Phys. Lett.} {\bf 89A}, 33 (1982)}
studied the three-dimensional 3-state Potts model with the usual ferromagnetic
nearest-neighbor interaction but also with an antiferromagnetic next-nearest-neighbor
interaction; they found a second order PT at $T_c = 0.845 \pm 0.005$.
This has recently been confirmed 
\ref\FMPRT{L. A. Fernandez, E. Marinari, G. Parisi, S. Roncolini and A. Tarancon,
``The 3d $Z_3$ Spin Model and the Deconfinement Transition in QCD: a Problem of
Universality'', ROM2F-88-032 Preprint (October 1988)}.
For larger q in three dimensions, there is the work by Ono and Ito
\ref\OnoIto{I. Ono and K. Ito, {\it J. Phys. C} {\bf 15}, 4417 (1982)}
showing that $T_c \approx 1.58$ for $q = 4$ and $T_c \approx 1.33$ for $q = 6$.
(They actually studied a lattice model for liquid crystals which is
completely equivalent to the Potts model for $q = 3,4$ or $6$.)

\bigskip
\line{\bf 2.2. Ising (q=2) model\hfill}

The earliest reference to a computer simulation of the Ising model appears to be that of
Fosdick \ref\Fosdick{L. D. Fosdick, {\it Bull. Am. Phys. Soc.} {\bf 2}, 239 (1957)}
who performed a Monte Carlo simulation of the two-dimensional
Ising model using $10^2, 20^2$ and $37^2$ lattices on the ILLIAC computer.
He found ``away from the Curie point, errors of less than 1\%'' with
the number of iterations varying ``up to 768 passes through the lattice'', and
observed that ``100 passes appeared to bring the system to `equilibrium' except in
the neighborhood of the Curie temperature where fluctuations cause erratic results''!
Thus began over thirty years of Ising computer simulations -- I shall just survey a few
highlights from the many papers produced.
A good history up to 1967 can be found in
\ref\Brush{S. G. Brush, {\it Rev. Mod. Phys.} {\bf 39}, 883 (1967)} and a thorough review
(of the 2d model)
in \ref\McCoyWu{B. M. McCoy and T. T. Wu, {\it The Two Dimensional Ising Model}
(Harvard U. P., Cambridge, Mass., 1973)}.
Landau, in 1976, produced comprehensive papers on finite size scaling for both the 2d
\ref\Landautwo{D. P. Landau, \PRB{13}, 2997 (1976)}
and 3d \ref\Landauthree{D. P. Landau, \PRB{14}, 255 (1976)} Ising models.
The end of the road of such ``brute force'' methods for obtaining critical
exponents was signaled when Monte Carlo Renormalization Group appeared on the
scene in 1979, and finally came in 1985 when a specially constructed processor for
performing Monte Carlo simulations on the 3d Ising model
\ref\SantaB{R. B. Pearson, J. L. Richardson and D. Toussaint, {\it J. Comput. Phys.}
{\bf 51}, 241 (1983);
M. N. Barber, R. B. Pearson, D. Toussaint and J. L. Richardson, \PRB{32}, 1720 (1985)}
failed due to a bad pseudo-random number generator
\ref\PRbad{G. Parisi and F. Rapuano, \PL{157}, 301 (1985)}
\ref\HCaBbad{A. Hoogland, A. Compagner and H. W. J. Bl\"ote, {\it Physica} {\bf 132A}, 593
(1985)}
\ref\BDaSbad{G. Bhanot, D. Duke and R. Salvador, \PRB{33}, 7841 (1986)}.
Interestingly, the second of these demonstrations of the failure of the Santa Barbara
Ising Model Processor (IMP) \HCaBbad\ made use of another special purpose machine -- the
Delft Ising System Processor (DISP)
\ref\HCaBbadtwo{A. Hoogland, J. Spaa, B. Selman and A. Compagner, {\it J. Comp. Phys.}
{\bf 51}, 250 (1983)}
, which was slower but had good random numbers!
Swendsen \ref\Swendsenone{R. H. Swendsen, \PRL{42}, 859 (1979)}
came up with a practical implementation of Ma's
\ref\Ma{S.-K. Ma, \PRL{37}, 461 (1976)}
suggestion for combining Monte Carlo simulations of statistical mechanical models with a
renormalization group
\ref\Wilson{K. G. Wilson, in {\it Recent Developments in Gauge Theories}, Cargese, 1979,
ed. G. t' Hooft (Plenum Press, New York, 1980)}
analysis of the critical properties, which he called Monte Carlo Renormalization Group
(MCRG).
In his original paper \Swendsenone, Swendsen looked at the 2d Ising model (and 3-state
Potts model), and then with Bl\"ote
\ref\Swendsentwo{H. W. J. Bl\"ote and R. H. Swendsen, \PRB{20}, 2077 (1979)}
investigated the 3d Ising model. These calculations were repeated five years later
\ref\Swendsenthree{R. H. Swendsen, \PRB{30}, 3866 (1984);
R. H. Swendsen, \PRB{30}, 3875 (1984)} culminating with the best numerical
work on the
3d Ising model so far
\ref\Edin{G. S. Pawley, R. H. Swendsen, D. J. Wallace and K. G. Wilson, \PRB{29}, 4030
(1984)} which determined the critical coupling to six decimal places:
$\beta_c = 0.221654(6)$, and the thermal and magnetic exponents to
three decimal places: $\nu = 0.629(4)$ and $\eta = 0.031(5)$.
Bl\"ote and Swendsen \ref\BaSI{H. W. J. Bl\"ote and R. H. Swendsen, \PRB{22}, 4481 (1980)}
also investigated the 4d Ising model with MCRG, confirming the critical coupling
$\beta_c = 0.29962$ and verifying the mean-field values of the critical exponents.

\bigskip
All these Monte Carlo values of $T_c$, for Potts and Ising models,
are summarized in Table 3.
\medskip
\vbox{
\bigskip
\input sptab3
\bigskip
}

\bigskip
\line{\bf 2.3. O(N) models\hfill}

Shenker and Tobochnik
\ref\SaT{S. H. Shenker and J. Tobochnik, \PRB{22}, 4462 (1980)}
performed the first computer simulation of the two-dimensional O(3) model,
using MCRG, showing that it does not have a phase transition.
Soon after, Martinelli {\it et al}
\ref\MPP{G. Martinelli, G. Parisi and R. Petronzio, \PL{100}, 485 (1981)}
found no phase transition using standard Monte Carlo methods.
Later Fox {\it et al}
\ref\FGMO{G. Fox, R. Gupta, O. Martin and S. Otto, \NP{205}, 188 (1982)}
estimated the mass gap for both O(2) and O(3) models.
Then Fukugita and Oyanagi \ref\FaO{M. Fukugita and Y. Oyanagi, \PL{123}, 71 (1983)}
calculated the mass gap and susceptibility for O(2), O(3), O(4) and O(5), verifying
the scaling behavior of those quantities with N.
There the situation would have rested had not Seiler {\it et al}
\ref\SSPL{E. Seiler, I. O. Stamatescu, A. Patrascioiu and V. Linke, \NP{305}, 623 (1988)}
claimed that their ``high precision numerical study'' of O(2) and O(4) gave a
second order phase transition for both. We shall discuss the case of O(2) (i.e. XY model)
below. This incorrect claim for O(4) was dismissed almost immediately by Heller
\ref\Hell{U. M. Heller, \PRL{60}, 2235 (1988)} using an improved MCRG method, and by
Bernreuther and G\"ockeler
\ref\BaG{W. Bernreuther and M. G\"ockeler, \PL{214}, 109 (1988)}
using Binder's method for locating critical points.

\bigskip
\line{\bf 2.4. XY (N=2) model\hfill}

The Kosterlitz and Thouless (KT) theory for the 2d XY model
\ref\KTK{J. M. Kosterlitz and D. J. Thouless, {\it J. Phys. C} {\bf 6}, 1181 (1973);
J. M. Kosterlitz, {\it J. Phys. C} {\bf 7}, 1046 (1974)}
predicts an essential singularity at a critical temperature $T_c = 1.35$ where the
correlation length and susceptibility diverge with exponents $\nu = {1 \over 2}$
and $\eta = {1 \over 4}$, respectively.
The first Monte Carlo study of the 2d XY model was by Tobochnik and Chester in
1979 \ref\TaC{J. Tobochnik and G. V. Chester, \PRB{20}, 3761 (1979)}, in which they
found a peak in the specific heat at $T = 1.02$, 15\% above a critical temperature
$T_c = 0.89$.
The phenomenon of vortex pairs unbinding, which is associated with
such a KT phase transition, was observed.
Unfortunately they could not obtain an unambiguous fit to their data and so ended up with
two fits for the exponent $\nu$ (extracted from the divergence of the susceptibility):
$0.70$ with $T_c = 0.89$; $0.48$ with $T_c = 0.92$.
They also obtained two probable values for $\eta$: $0.254$ at $T = 0.90$;
$0.243$ at $T = 0.85$.
Later, with finite size scaling, Fernandez {\it et al}
\ref\Fer{J. F. Fernandez, M. F. Ferreira and J. Stankiewicz, \PRB{34}, 292 (1986)}
managed to obtain $T_c = 0.89 \pm 0.02$, $\nu = 0.5 \pm 0.1$ and $\eta = 0.24 \pm 0.03$.
Next, of course, Seiler {\it et al} \SSPL\ claimed that the XY model phase transition
was actually second order because they obtained a better fit for this compared to their
KT fit. This KT fit gave $T_c = 0.7463 \pm 0.0002$, $\nu = 1.689 \pm 0.004$ and
$\eta = 0.47 \pm 0.19$.
They were quickly shown to be wrong by Gupta {\it et al}
\ref\XY{R. Gupta, J. DeLapp, G. G. Batrouni, G. C. Fox, C. F. Baillie and J. Apostolakis,
\PRL{61}, 1996 (1988)} with a high statistics, large lattice study, and by
Bernreuther and G\"ockeler \BaG\ using Binder's method for locating critical points.
In fact, Gupta {\it et al} did more than just prove Seiler {\it et al} wrong -- they
cleared up the situation by finding $four$ stable fits to their XY data:
three spurious minima and the true minimum!
The basin of attraction of the first spurious minimum
($T_c = 0.701$, $\nu = 1.460$ and $\eta = 0.53$)
includes the results of \SSPL, while the second and third spurious minima
($T_c = 0.867$, $\nu = 0.701$ and $\eta = 0.32$;
$T_c = 0.911$, $\nu = 0.472$ and $\eta = 0.34$)
are consistent with the fits of \TaC.
The true minimum is $T_c = 0.8987(1)$, $\nu = 0.4995(11)$ and $\eta = 0.34$
(from the susceptibility; from the correlation length Gupta {\it et al} get
$T_c = 0.8967(2)$, $\nu = 0.5007(3)$ and $\eta = 0.34$).

The 3d XY model has also been analyzed -- Swendsen
\ref\SXY{R. H. Swendsen, \PRB{27}, 391 (1983)}, using MCRG, confirmed the critical
coupling (from series expansions) $\beta_c = 0.45393$.

\bigskip
\line{\bf 3. How to beat critical slowing down\hfill}

In order to beat CSD we must first ask: what causes critical slowing down?
The key feature about standard (Metropolis-like) Monte Carlo algorithms is that
the updates are {\it local}, that is, one spin at a time is updated.
Thus in a single step of the algorithm, ``information'' about the state of a spin
is transmitted only to its nearest neighbors. Now, in order for the system to reach a new
{\it effectively independent} configuration, this information must travel a distance of
order the (static or spatial) correlation length $\xi$.
As the information executes a random walk around the lattice, one would suppose that
the autocorrelation time $\tau \sim \xi^2$ near criticality.
This can be shown analytically to be correct for a free-field (Gaussian) model.
However, in general,
$$\tau \sim \xi^z, \eqno\eq$$
where $z$ is called the dynamical critical exponent.
Almost all numerical computer simulations of spin models
have measured $z \approx 2$ for local update algorithms (with one exception -- the
over-relaxed algorithm -- described in Sec. 4), as
we shall review in the next subsection before going on
to explain the new algorithms in Sec. 5 which reduce $z$ by
performing {\it nonlocal} spin updates.

Firstly, we should briefly explain how the standard (single-site)
Metropolis \Metro\ algorithm works
for completeness. In the Metropolis algorithm a new configuration
is generated by updating a single variable in the old configuration $\phi$ and
calculating the change in energy
$$\Delta S = S(\phi') - S(\phi).\eqno\eq$$
If $\Delta S \le 0$, the change is accepted; if $\Delta S > 0$ the change
is accepted with probability $\exp(- \Delta S)$.
In practice, on a computer, this is done by generating a pseudo-random number $r$ in the
interval [0,1) with uniform probability distribution and accepting the change
if $r < \exp(- \Delta S)$.
If the variable is updated by a large amount then the change in energy will be large
and the probability of the update being accepted will be small.
Conversely, if the variable is changed by a small amount then the changes will usually
be accepted but their small size will make exploration of phase space rather slow.
The lore is for a compromise: tuning the update size to get an acceptance rate of about
50\%.

\bigskip
\line{\bf 3.1. Measurements of $z$ using local update algorithms\hfill}

A great deal of effort has been spent verifying that $z \approx 2$ using local
update algorithms for a variety of spin models.

\bigskip
\line{\bf 3.1.1. 2d Potts\hfill}

Binder \Binder, as well as looking at the static model, investigated the dynamic 2d Potts
model for $q = 3,4,5,6$, finding $z = 1.9$ independent of $q$ within his accuracy.
The first computer calculation using the dynamical MCRG method was by
Tobochnik and Jayaprakash
\ref\TobJay{J. Tobochnik and C. Jayaprakash, \PRB{25}, 4893 (1982)} for
$q = 3$ to obtain $z = 2.7 \pm 0.4$.
Aydin and Yalabik
\ref\AY{M. Aydin and M. C. Yalabik, {\it J. Phys. A} {\bf 17}, 2531 (1984)}
repeated this calculation with a modified dynamical MCRG method and got the same
answer; then they went on to redo the calculation
\ref\AYY{M. Aydin and M. C. Yalabik, {\it J. Phys. A} {\bf 18}, 1741 (1985)}
using finite size scaling and getting $z = 2.2 \pm 0.1$.
Next de Arcangelis and Jan
\ref\Lucilla{L. de Arcangelis and N. Jan, {\it J. Phys. A} {\bf 19}, L1179 (1986)}
found $z = 2.43 \pm 0.15$ for $q = 3$, and $z = 2.36 \pm 0.20$ for $q = 4$,
using the modified dynamical MCRG method.
Then Jain \ref\Jain{S. Jain, {\it J. Phys. A} {\bf 20}, L607 (1987)} obtained
$z = 2.94 \pm 0.29$ for $q = 4$ with finite size scaling.
Finally Tang and Landau
\ref\TangLand{S. Tang and D. P. Landau, \PRB{36}, 567 (1987)} determined
$z = 2.17 \pm 0.04$ for $q = 2,3$ and $4$ with extensive runs using finite size scaling.
These values of $z$ are summarized in Table 4.
\medskip
\vbox{
\bigskip
\input sptab4
\bigskip
}
No one appears to have done any numerical computer simulations in order to calculate
$z$ for the Potts model in more than two dimensions.

\bigskip
\line{\bf 3.1.2. 2d Ising\hfill}

The earliest computer work
\ref\SSB{T. Schneider, E. Stoll and K. Binder, \PRL{29}, 1080 (1972);
E. Stoll, K. Binder and T. Schneider, \PRB{8}, 3266 (1973)}
on $z$ for the 2d Ising model
did not use finite size scaling and got the wrong answer $1.85 \pm 0.10$.
The first reliable computer calculation was by Tobochnik {\it et al}
\ref\TSaC{J. Tobochnik, S. Sarker and R. Cordery, \PRL{46}, 1417 (1981)}
in which they devised the dynamical MCRG method and used it to measure
$z = 2.22 \pm 0.13$.
This was repeated by Katz {\it et al}
\ref\KGaL{S. L. Katz, J. D. Gunton and C. P. Liu, \PRB{25}, 6008 (1982)}
giving $z = 2.23$. At the same time d'Auriac {\it et al}
\ref\AMaR{J. C. Angles d'Auriac, R. Maynard and R. Rammal, {\it J. Stat. Phys.}
{\bf 28}, 307 (1892)}
used finite size scaling with a clever update algorithm to obtain $z = 2.0 \pm 0.1$.
Later Jan {\it et al}
\ref\JMaS{N. Jan, L. L. Moseley and D. Stauffer, {\it J. Stat. Phys.} {\bf 33}, 1
(1983)} came up with the modified dynamical MCRG method and obtained $z = 2.12 \pm 0.06$
with it (they also looked at 3d -- see below).
With more computer time Kalle \ref\Kalle{C. Kalle, {\it J. Phys. A} {\bf 17}, L801 (1984)}
repeated this and obtained $z = 2.14 \pm 0.02$ (as well as the value for 3d -- see below).
Another careful, extensive study by Williams
\ref\Will{J. K. Williams, {\it J. Phys. A} {\bf 18}, 49 (1985)}
yielded $z = 2.13 \pm 0.03$.
Finally, using a more traditional finite size scaling technique Miyashita and Takano
\ref\MaT{S. Miyashita and H. Takano, {\it Prog. Theor. Phys.} {\bf 73}, 1122 (1985)}
obtained $z = 2.2$.

\bigskip
\line{\bf 3.1.3. 3d Ising\hfill}

Early calculations of $z$ for the 3d Ising model were by Chakrabarti {\it et al}
\ref\CBaS{B. K. Chakrabarti, H. G. Baumg\"artel and D. Stauffer, {\it Z. Phys. B} {\bf 44},
333 (1981)} using finite size scaling to get $2.17 \pm 0.06$,
and by Yalabik and Gunton
\ref\YaG{M. C. Yalabik and J. D. Gunton, \PRB{25}, 534 (1982)} using dynamical
MCRG to obtain $z = 2.08$.
Later Jan {\it et al}
\JMaS\ measured $z = 2.11 \pm 0.03$, and Kalle \Kalle\ found $z = 1.965 \pm 0.010$.
The Santa Barbara group
\ref\SantaBB{R. B. Pearson, J. L. Richardson and D. Toussaint, \PRB{31}, 4472 (1985)}
also measured $z$ -- on their Ising Model Processor -- obtaining $1.99 \pm 0.02$ using
finite size scaling, but this result cannot be trusted in the light of later revelations
about their random numbers (see above).
Finally, Wansleben and Landau
\ref\WaL{S. Wansleben and D. P. Landau, {\it J. Appl. Phys.} {\bf 61}, 3968 (1987)}
obtained $z = 2.03 \pm 0.04$ using finite size scaling, and argued that, since $z$ should be
greater than $2$ (on theoretical grounds from $\epsilon$ expansions),
\Kalle\ has underestimated it and \SantaBB\ $is$ wrong.

\bigskip
These values of $z$, for both the 2d and 3d Ising models, are summarized in Table 5.
\medskip
\vbox{
\bigskip
\input sptab5
\bigskip
}
\bigskip

\bigskip
\line{\bf 3.1.4. 2d XY\hfill}

Fernandez {\it et al} \Fer\ measured $z = 2$ during their 2d XY model simulation.

\bigskip
\line{\bf 4. Over-relaxed algorithm\hfill}

The over-relaxed algorithm was invented by Adler
\ref\Adler{S. L. Adler, \PRD{23}, 2901 (1981); S. L. Adler, \PRD{37}, 458 (1988)},
and first applied to lattice gauge theories by Creutz
\ref\Creutz{M. Creutz, \PRD{36}, 515 (1987)}
and by Brown and Woch
\ref\BaW{F. R. Brown and T. J. Woch, \PRL{58}, 2394 (1987)}.
Adler showed analytically that this algorithm reduces the dynamical
critical exponent from $2$ to $1$ for the Gaussian model.
The over-relaxed algorithm is applicable to any theory where the variables are
elements of a group larger than $Z_2$ (Ising).
It is a modification of the single-site Metropolis \Metro\ Monte Carlo algorithm.
The basic idea is to pick a trial change for a given variable in a region of phase
space which is as far as possible from the old value but which leaves the energy
the same. Then the probability of acceptance is unity, yet the subsequent
configurations decorrelate much faster than in the Metropolis algorithm.
However, the fact that the energy remains constant means that the algorithm is
deterministic and may not explore the whole phase space. Therefore, to
ensure ergodicity, a standard Metropolis sweep of the lattice is done after
every few over-relaxed sweeps.
Brown and Woch give a rather nice picture in order to explain how this over-relaxed
algorithm fits into the scheme of things; we now present a slightly different
version of this.
Given a quadratic action (for simplicity)
$$S(x) = a (x - x_0)^2 + c, \eqno\eq$$
the following over-relaxed update satisfies detailed balance with respect to $S$:
$$x' = \omega x_0 + (1-\omega) x + w(2-\omega) a^{-1/2} \eta, \eqno\eq$$
where $\eta$ is unit Gaussian noise, $\omega$ is the over-relaxation parameter
and $x_0$ is the environment (i.e. sum of nearest neighbor $x$'s) of $x$.
(The analogy is with improving the convergence of iterative linear algebra solvers
via the introduction of an over-relaxation parameter, that is, changing an
equation of the form
$$x^{(n)} = M r^{(n-1)} + x^{(n-1)} \eqno\eq$$
into
$$x^{(n)} = \omega M r^{(n-1)} + (1-\omega) x^{(n-1)} \eqno\eq$$
.)
We see that the resulting algorithm depends on $\omega$ as depicted in Figure 1.
For $\omega = 0$ the system does not change at all, for $0 < \omega < 1$
(under-relaxation) we get the Metropolis algorithm (with $\omega$ corresponding to the
update size), for $\omega = 1$ we get the conventional heat bath algorithm, for
$1 < \omega < 2$ (over-relaxation) we get the general over-relaxed algorithm,
and for $\omega = 2$ we get the specific microcanonical over-relaxed algorithm
(with deterministic constant-energy evolution) which we want.
(The algorithm diverges if $\omega \not\in [0,2]$.)
The over-relaxed algorithm was invented just in time to be used for one
state-of-the-art spin model computer simulation (see next subsection),
before being overtaken by the so-called cluster algorithms.

\bigskip
\line{\bf 4.1. Measurements of $z$ using over-relaxed algorithm\hfill}

The only spin model calculation in which the over-relaxed algorithm was used
was that of Gupta {\it et al} \XY. They used it in conjunction with the standard
Metropolis algorithm. With 8 sweeps of over-relaxed and 2 sweeps of Metropolis they
obtained $z = 1.48$, and with 15 OR sweeps and 2 Met sweeps they got $z = 1.20$.
Thus the over-relaxed/Metropolis algorithm is indeed the best {\it local} update
algorithm - allowing one to reduce the dynamical critical exponent from $2$ to
close to $1$.
However, one can do better with {\it nonlocal} cluster algorithms.

\bigskip
\line{\bf 5. Cluster algorithms\hfill}

Let us review the reasoning of why nonlocal updating should ameliorate CSD.
The physics of critical slowing down tells us that the slow modes in the theory or model
are the long-wavelength modes, if the updating is purely local.
The natural physical solution is therefore to speed up these modes by
{\it collective-mode} (nonlocal) updating. Of course the catch is that one must
decide exactly which collective modes to speed up and how. Moreover the algorithm should
not be so nonlocal that its computational cost outweighs the reduction in CSD.
There are at least three ways to implement nonlocal, collective-mode updating:
\medskip
\item{i.} Fourier acceleration
\smallskip
\item{ii.} Multigrid Monte Carlo
\smallskip
\item{iii.} Cluster algorithms
\medskip
\noindent
We shall not discuss Fourier acceleration or Multigrid Monte Carlo here.
In any case, the cluster algorithms are easier to program and seem to
give as good performance.

The aim of the cluster algorithms is to find a suitable collection of spins which
can be flipped with relatively little cost in energy.
Note that we could obtain nonlocal updating very simply by using the standard Metropolis
Monte Carlo algorithm to flip randomly selected bunches of spins. But then the
acceptance would be tiny. Therefore we need a method which picks sensible bunches or
clusters of spins to be updated. Currently there are three such methods available --
Swendsen and Wang, Niedermayer, and Wolff.

All these algorithms have in common the problem of identifying clusters of spins --
we explain the best known sequential algorithm for this (by Hoshen and Kopelman)
and give two parallel algorithms.
We are not looking, here, at the measurement part of lattice spin model computer simulations
but we should note in passing that the clusters built by cluster algorithms can be used
for measuring physical observables and in fact yield values with smaller statistical errors
than are obtained from the spins directly.

\bigskip
\line{\bf 5.1. Swendsen and Wang\hfill}

Cluster algorithms became important for spin models when Swendsen and Wang (SW)
published their algorithm for Potts models
\ref\SW{R. H. Swendsen and J.-S. Wang, \PRL{58}, 86 (1987)}.
This is based on an equivalence between configurations of Potts spins and
configurations of clusters in a percolation problem, which was discovered by
Kasteleyn and Fortuin some time ago
\ref\KF{P. W. Kasteleyn and C. M. Fortuin, {\it J. Phys. Soc. Jpn. Suppl.} {\bf 26s},
11 (1969); C. M. Fortuin and P. W. Kasteleyn, {\it Physica} {\bf 57}, 536 (1972)}.
Sweeny \ref\Sweeny{M. Sweeny, \PRB{27}, 4445 (1983)} also noticed
Kasteleyn and Fortuin's work and came up with a different algorithm for simulating the
q-state Potts model for arbitrary (even non-integer) q. He achieved this by
simulating the equivalent percolation clusters. However, Sweeny's algorithm is complicated
to program and is restricted to two dimensions, whereas SW is relatively easy to program and
works in any dimension (but only for integer q). Moreover, SW has been quickly taken up and
generalized to continuous spin models by Niedermayer and Wolff (see below).
We shall firstly describe SW and summarize applications of it, then outline
Niedermayer, and finally concentrate on Wolff as it shows the most promise.

To explain SW, we begin with the partition function for the Potts model
$$Z = Tr_{\{\sigma\}} e^{\beta \sum_{<ij>} (\delta_{\sigma_i,\sigma_j} - 1)}. \eqn\EqSWo$$
We map it onto a general percolation model {\it \`a la} Kasteleyn and Fortuin, as
follows.
Firstly remove the interaction between sites $l$ and $m$
$$S_{l,m} = \beta \sum_{<i,j> \not= <l,m>} (\delta_{\sigma_i,\sigma_j} - 1). \eqno\eq$$
Then define the restricted sums
$$Z_{l,m}^{same} = Tr_{\{\sigma\}} e^{S_{l,m}} \delta_{\sigma_l,\sigma_m} \eqno\eq$$
and
$$Z_{l,m}^{diff} = Tr_{\{\sigma\}} e^{S_{l,m}} (1 - \delta_{\sigma_l,\sigma_m}) \eqno\eq$$
so that
$$Z = Z_{l,m}^{same} + e^{-\beta} Z_{l,m}^{diff}. \eqn\EqSWt$$
Now introduce
$$Z_{l,m}^{indep} = Tr_{\{\sigma\}} e^{S_{l,m}} = Z_{l,m}^{same} + Z_{l,m}^{diff}. \eqno\eq$$
and rewrite \EqSWt\ as
$$Z = (1 - e^{-\beta}) Z_{l,m}^{same} + e^{-\beta} Z_{l,m}^{indep}. \eqno\eq$$
Since the first term contains the restriction that the spins on sites $l$ and $m$ are
the same, and the second term places no restrictions on the spins, the weighting
factors can be interpreted in terms of the probability $p = 1 - e^{-\beta}$ of a
bond between the two sites.
This procedure is repeated for all interactions in \EqSWo\ giving a configuration of
bonds, from which clusters are formed. For each cluster, all the spins are in one of the
q states independently of the states of the other clusters. If there are $N_c$ clusters,
this gives Kasteleyn and Fortuin's result
$$Z = Tr_{\{bonds\}} p^b (1-p)^n q^{N_c}, \eqno\eq$$
where $b$ is the number of bonds and $n$ is the number of interactions that do not form
a bond. Thus we have represented the spin system as a configuration of bonds forming
clusters. Notice that some of the original clusters of spins may be split into
smaller clusters, since bonds occur only with probability $p = 1 - e^{-\beta}$ between
sites with the same spin.
SW's algorithm consists of assigning a new random Potts spin value to each cluster
independently, with all the sites in the cluster getting this value.
By erasing bonds we are left with a new Potts spin configuration. This new
configuration can differ substantially from the original one since large clusters
can be changed in a single step. Thus we have achieved a nonlocal move which is what is
required to reduce CSD. SW \SW\ demonstrate that their algorithm satisfies detailed
balance.

\bigskip
\line{\bf 5.1.2. Applications -- measurements of $z$\hfill}

In their original paper \SW, SW measured the dynamical critical exponent for the 2d and 3d
Ising model, obtaining $z = 0.35(1)$ and $z = 0.75(1)$ respectively, as well as for
the 2d 3-state Potts model, getting $z = 0.6(1)$.
This dramatic result means that the SW algorithm is over 1000 times more
efficient than standard Metropolis (for which $z \approx 2$) on a $512^2$ lattice.
Montvay {\it et al}
\ref\Mont{I. Montvay, G. M\"unster and U. Wolff, \NP{305}, 143 (1988)}
used SW to study the 4d Ising model. On $12^4$ and $24^4$ lattices, where the (spatial)
correlation lengths are around $2$ and $4$, they found that the relative errors
using SW were down by factors of $2-2.5$ and $4-5$, respectively, over standard Metropolis.
This means that the autocorrelation lengths are down by roughly $5$ and $20$,
respectively, so that the total gain in computational speed is around an order of
magnitude.
Despite this extraordinary performance on Ising/Potts models, the SW algorithm is not
very good for O(N) models. In fact, for the 2d XY model it only reduces CSD on the low
temperature side of criticality; on the high temperature side there is still 
significant CSD
\ref\Sokal{A. D. Sokal, in {\it Springer Proceedings in Physics} {\bf 33}, 6
(Springer-Verlag, Berlin, 1988)}.
Therefore we turn to generalizations of the SW algorithm,
specifically aimed at O(N) models.

\bigskip
\line{\bf 5.2. Niedermayer\hfill}

Niedermayer \ref\N{F. Niedermayer, \PRL{61}, 2026 (1988)} generalized SW in the following
way. Start with the general action
$$S = \beta \sum_{<i,j>} E_l(\sigma_i,\sigma_j), \eqno\eq$$
where $E_l$ is the term in the total energy corresponding to a given pair of sites
on which live the continuous spins $\sigma_i$ and $\sigma_j$.
Assume that $S$ is invariant under a global symmetry transformation $g$,
$$\sigma_i \rightarrow \sigma_i' = g \sigma_i:
\quad E_l(\sigma_i,\sigma_j) = E_l(\sigma_i',\sigma_j'). \eqn\EqNo$$
Now introduce bonds on links $l$ joining sites $i$ and $j$ with probability
$p_l(E_l(\sigma_i,\sigma_j))$, where the probability function $p_l$ is
arbitrary. Sites connected by bonds form clusters as in SW.
The probability of obtaining a particular configuration of clusters $C$ from a
configuration of spins $S$ is
$$w(S \rightarrow C) = \sum_{B(C)} \prod_{l \in B} p_l(E_l) \prod_{l \not\in B} (1-p_l(E_l)),
\eqno\eq$$
where products are taken over links with and without bonds $B(C)$ respectively.
We then imagine clusters as separate subsystems and make an independent global
symmetry transformation on each one. Since bond probabilities are arbitrary we must
calculate the change in energy for transforming each cluster and accept/reject just
as in standard (single-site) Metropolis.
As $E_l = E_l'$ when both end points of the link belong to the same cluster \EqNo,
the change in energy $\Delta E$ comes only from the boundary $\partial C$ of the cluster
and is given by
$$\Delta E = \prod_{\partial C} e^{-\beta(E_l' - E_l)}
{(1-p_l(E_l')) \over (1-p_l(E_l))}. \eqn\EqNt$$
Niedermayer \N\ shows that this satisfies detailed balance.

In order to compare Niedermayer's algorithm with SW, let us make a special choice of the
bond probabilities:
$$p_l(E_l) = \cases{ 1 - e^{\beta(E_l-E_0)}, &$E_l < E_0$\cr
                     0, &$E_l \geq E_0$\cr}$$
with $E_0$ a free parameter.
Since $p_l(E_l)$ decreases as $E_l$ increases (up to $E_0$), the smaller the energy of
the link, the larger the probability that the two corresponding sites belong to the
same cluster. Thus these clusters should be similar to the real clusters in the
physical system. For $E_0 \leq E_{min}$, $p_l = 0$ and each cluster consists of only
a single site -- this is just the standard Metropolis local updating procedure and
\EqNt\ reduces to $\Delta E = e^{-\beta(E_l' - E_l)}$ for each link.
With increasing $E_0$ the average size of the clusters grows.
For $E_0 \geq E_{max}$, $p_l = 1 - e^{\beta(E_l - E_0)}$ and all factors cancel in
\EqNt\ which means that there is no interaction between clusters and one can choose
the new configuration at random -- thus we have recovered SW.
In fact, for the Ising model the choice $E_0 = J = 1 = E_{max}$ gives
$p = 1 - e^{-2\beta}$ for parallel spins and $p = 0$ for antiparallel ones,
which are the probabilities used in SW. (The extra factor of two here is due to
the action for Ising \EqIsing\ involving the product of the spins whereas the
Potts action contains a delta-function \EqSWo.)
The advantage of Niedermayer's algorithm is that one can tune $E_{min} < E_0 < E_{max}$
to get the best average cluster size for minimizing the relaxation time.

\bigskip
\line{\bf 5.2.2. Applications -- measurements of $z$\hfill}

Niedermayer \N\ tested his algorithm on the 2d XY model. He didn't gather enough data
to determine $z$ but he did see a ``considerable decrease in relaxation time''.
The disadvantage with his algorithm is, of course, the fact that the change in energy
must be calculated around the boundary of each cluster, which is tricky to program and
time consuming. Fortunately there is another generalization of SW which avoids this
difficulty.

\bigskip
\line{\bf 5.3. Wolff\hfill}

Wolff \ref\W{U. Wolff, \PRL{62}, 361 (1989)} generalized SW for O(N) models by inventing
an appropriate spin flip operation. Consider the O(N) model partition function
$$Z = Tr_{\{{\bf\sigma}\}} e^{\beta \sum_{<ij>} {\bf\sigma_i . \sigma_j}}, \eqno\eq$$
where ${\bf\sigma_i}$ are unit vectors in $\Re^N$.
Generalize the spin flip operation $\sigma_i \rightarrow -\sigma_i$ for the Ising model
($N=1$) in the following way. For $N \geq 2$ and any ${\bf r} \in S_{N-1}$ define
the generalized spin flip as the reflection of ${\bf\sigma_i}$ with respect to the
hyperplane orthogonal to ${\bf r}$:
$${\bf R}({\bf r}){\bf\sigma_i} = {\bf\sigma_i} - 2 ({\bf\sigma_i . r}){\bf r}. \eqn\EqWo$$
Clearly, the action is invariant under this global ${\bf R}$ transformation since
$$[{\bf R}({\bf r}){\bf\sigma_i}] . [{\bf R}({\bf r}){\bf\sigma_j}] =
{\bf\sigma_i . \sigma_j}. \eqno\eq$$
An elementary cluster update consists of the following steps
\medskip
\item{i.} choose random reflection ${\bf r} \in S_{N-1}$ and random site $i$ as first
point of cluster $c$ to be built
\smallskip
\item{ii.} flip ${\bf\sigma_i} \rightarrow {\bf R}({\bf\sigma_i})$ and mark $i$
\smallskip
\item{iii.} visit all links $l$ connecting $i \in c$ to its nearest neighbors $j$ and
activate bond $l = <ij>$ with probability
$$\eqalign{p({\bf\sigma_i,\sigma_j})
&=\ 1 - e^{min[0,\beta{\bf\sigma_i}.({\bf 1 - R}({\bf r})){\bf\sigma_j}]}\cr
&=\ 1 - e^{min[0,-2\beta({\bf r.\sigma_i})({\bf r.\sigma_j})]}\cr
}\eqn\EqWt$$
and if this happens ${\bf \sigma_j}$ is flipped, and $j$ is marked and joined to $c$
\item{iii.} continue iteratively for all bonds leading to unmarked neighbors of
newly adjoined sites until process stops
\medskip
\noindent
Wolff \W\ shows that this satisfies detailed balance.

For the Ising case the bond activation probabilities are the same as in SW:
$p = 1 - e^{-2\beta}$ for parallel spins and $p = 0$ for antiparallel ones.
The connection between the two algorithms can be made clearer by consideration
of a less efficient but equivalent implementation of Wolff's algorithm.
This consists of first activating all bonds of the lattice with the
appropriate probabilities and decomposing all the sites into clusters. These
are the SW clusters. Then flip each cluster independently with probability
${1 \over 2}$, rather than selecting one cluster and flipping it alone.

The adequacy of discrete flips for continuous spin models is ensured because the
axis of reflection is chosen randomly.
If we choose ${\bf r}$ and build a cluster in a region where spins are nearly
orthogonal to ${\bf r}$ then \EqWo\ is a small deformation and bond probabilities
\EqWt\ typically lead to small clusters. On the other hand, if ${\bf r}$ is
parallel or antiparallel to the local spin direction then the cluster tends to
grow until it reaches an energetically acceptable surface -- thus mimicking
real clusters.

Wolff's algorithm is better than Niedermayer's because it doesn't involve the
calculation of the change in energy on the boundary of a cluster in order to flip it --
instead this is built into the bond probabilities
(by restricting global symmetry transformations to reflections in a hyperplane)
as it is in SW. In fact, Wolff complains that when he tried Niedermayer's algorithm
he found that ``either the whole lattice ends in one cluster, or large
nontrivial clusters rarely move because of their surface energy'' \W.
Presumably this could be ameliorated by tuning the parameter $E_0$ in
Niedermayer's algorithm. In any case, Wolff's algorithm is easier to
program and gives as good if not better results than both Niedermayer's and SW --
therefore it is the best currently available algorithm for efficiently
updating (discrete and continuous) spin models.

\bigskip
\line{\bf 5.3.2. Applications -- measurements of $z$\hfill}

Wolff \W\ tested his algorithm on the 2d O(N) models with $N = 1,2,3$.
For $N=1$ (Ising) and $N=2$ (XY) he measured the autocorrelation length
on lattices of size $32^2, 64^2, 128^2$.
For the Ising he obtained $2.3(3), 2.3(3), 2.7(5)$ and for the XY
$3.6(7), 2.4(6), 2.1(6)$ -- both of which are consistent with $z \approx 0$ !
For the O(3) model which does not have a phase transition, Wolff ran at
a temperature which leads to a (spatial) correlation length of around 20
and found that the distribution of the number of spins in a cluster has
a cut-off at the expected maximum cluster size. Thus the Wolff clusters do
indeed mimic the real ones, and moreover his algorithm eliminates CSD altogether.
Recently Wolff has used his algorithm for detailed studies of both
the XY \ref\WXY{U. Wolff, ``Collective Monte Carlo Updating in a High Precision Study
of the XY Model'', DESY 88-176 Preprint (December 1988)}
and O(3) \ref\WO{U. Wolff, ``Asymptotic Freedom and Mass Generation in the O(3)
Nonlinear Sigma Model'', DESY 89-021 Preprint (February 1989)} models.
For the XY model, no CSD is seen with correlation lengths up to $89$.
More importantly, from the physics point of view, Wolff finds a Kosterlitz-Thouless
transition, not a second order one as was erroneously claimed (see Sec. 2.4 above).
Similarly, for the O(3) model the autocorrelation length is less than one.

\bigskip
\line{\bf 6. Cluster identification\hfill}

All cluster algorithms have in common the problem of identifying the clusters of spins.
The best known sequential algorithm for this is by Hoshen and Kopelman.
An effective parallel algorithm is the so-called ``ants in the labyrinth'' algorithm.
However, this identifies the clusters one at a time so we present our parallel algorithm,
``self-labeling'', which does them all at once.

\bigskip
\line{\bf 6.1. Hoshen and Kopelman labeling\hfill}

Hoshen and Kopelman \ref\HK{J. Hoshen and R. Kopelman, \PRB{14}, 3438 (1976)}
invented a sequential cluster identification algorithm which gives each cluster a unique
label and counts the number of sites it contains. Each site $i$ belonging to
cluster $\alpha$ is assigned a cluster label $m_t^\alpha$.
The cluster $\alpha$ may be assigned several cluster labels, which are given
as a set of natural numbers
$$\{ m_1^\alpha, m_2^\alpha, ..., m_s^\alpha, ..., m_t^\alpha, ...\}. \eqn\HKo$$
In this set only one number is regarded as a {\it proper} cluster label
which we designate as $m_s^\alpha$. This is the smallest number of the set in \HKo.
The following set of integers provides the connections between the $m_t^\alpha$
labels:
$$\{ N(m_1^\alpha), N(m_2^\alpha), ..., N(m_s^\alpha), ..., N(m_t^\alpha), ...\}. \eqn\HKt$$
Here only $N(m_s^\alpha)$ is positive and denotes the number of sites in the cluster.
The remaining members of the set in \HKt\ are negative integers, providing the links between
the other $m_t^\alpha$ labels and the proper label $m_s^\alpha$.
To be specific, if a site with label $m_t^\alpha$ is actually next to the
site which has the proper cluster label $m_s^\alpha$ then
$$N(m_t^\alpha) = -m_s^\alpha. \eqn\Ho$$
However, if a site with label $m_p^\alpha$ is not directly joined to the
site with the proper cluster label $m_s^\alpha$ but is joined to the site
with label $m_t^\alpha$ then
$$N(m_p^\alpha) = -m_t^\alpha. \eqn\Ht$$
Thus to get the proper cluster label for this site we have to through
$m_t^\alpha$ using \Ht\ then \Ho.
Similarly for higher-order indirections - fortunately in most cases this hierarchy extends
to one or two levels only.

This works in practice as follows. We sweep through sites of the lattice looking at
neighbors in the negative directions
(there are $d$ of them for a $d$-dimensional lattice).
If none of the neighbors of site $i$ have been labeled then $i$ is assigned a new label
$S_i = m_t^\alpha$, and $N(m_t^\alpha) = 1$.
If there is only one labeled neighbor, at site $n$ with label $S_n = m_r^\alpha$, say, then
$i$ gets the same label as $n$: $S_i = m_r^\alpha$, and $N(m_t^\alpha) = N(m_t^\alpha)+1$.
The good feature of this cluster labeling technique becomes apparent when the site $i$
links two or more previously labeled cluster fragments into a single cluster, that is,
when two or more of its neighbors are labeled.
No site belonging to any of these cluster fragments is relabeled (so that once a site is
labeled it retains this label throughout the labeling process) -- instead the readjustments
occur within the $N(m_t^\alpha)$ sets.
The number of readjusted $N(m_t^\alpha)$'s for site is equal to the number of
coalescing cluster fragments at that site. Let us assume that the neighboring sites
belong to clusters $\alpha, \beta, \gamma, ...$ which have proper cluster labels
$m_s^\alpha, m_s^\beta, m_s^\gamma, ...$ with $m_s^\alpha$ being the smallest.
These clusters coalesce at $i$ to form a single larger cluster so we set
label $S_i = m_s^\alpha$ and readjust
$$\eqalign{N(m_s^\alpha) &=\ N(m_s^\alpha) + N(m_s^\beta) + N(m_s^\gamma) + ... + 1\cr
N(m_s^\beta) &=\ - m_s^\alpha\cr
N(m_s^\gamma) &=\ - m_s^\alpha\cr
...\cr
}\eqno\eq$$
The only tricky part is finding the proper cluster labels for the neighboring clusters.
This is done in the following way: we start with the cluster label of the neighboring
site $S_n$ setting
{\settabs 10 \columns
$$\eqalign{r &=\ S_n\cr
t &=\ -N(r)\cr}$$
\+if $t < 0$ then $r$ is the proper cluster label so stop, else \cr
\+&set \cr
$$\eqalign{r &=\ t\cr
t &=\ -N(r)\cr}$$
\+&if $t < 0$ then $r$ is the proper cluster label so stop, else \cr
\+&&repeat setting \cr
$$\eqalign{r &=\ t\cr
t &=\ -N(r)\cr}$$
\+&&until $t < 0$ and $r$ is the proper cluster label and finally set \cr
$$N(S_n) = -r$$
\+&&which eliminates high-order levels in hierarchy of indirections. \cr
}
\smallskip
\noindent
As a final stage in the algorithm we can pass through the lattice a second time
copying each proper cluster label $m_s^\alpha$ to all sites belonging to that cluster.
This makes it easier to pick out the clusters for updating.

\bigskip
\line{\bf 6.2. Parallel labeling\hfill}

There are several parallel algorithms for cluster labeling, we just describe two of them.

\bigskip
\line{\bf 6.2.1. Ants in the labyrinth\hfill}

The best known parallelizable cluster identification algorithm is ``ants in the
labyrinth'' which is well described by Dewar and Harris
\ref\Ants{R. Dewar and C. K. Harris, {\it J. Phys. A} {\bf 20}, 985 (1987)}.
An ant is put somewhere on the lattice and notes what spin that site has.
At the next time-step this ant places children on each of the neighboring
sites which have the same spin. The children then proceed to reproduce likewise
until the entire cluster of such spins is populated. The parallel feature is that at
each time-step the whole of the frontier of the ant population on the cluster advances
at once. The ant population at any time is stored as a logical array (TRUE = ant present;
FALSE = no ant present), so this algorithm executes very quickly on bit-serial SIMD
computers like the AMT DAP and TMC Connection Machine.
The cluster is thus identified as a set of logicals (the ants)
whose size is the number of sites in the cluster. Once the cluster is identified it is
removed from the configuration. We repeat this ant-colony birth, growth and death
cycle until all the clusters have been identified and removed.
Of course, this algorithm is not fully parallel, because clusters are identified one by
one, so we turn to our algorithm, which is.

\bigskip
\line{\bf 6.2.2. Self-labeling\hfill}

The following algorithm is the one I immediately thought of when faced with the
problem of identifying clusters in parallel. It is so simple that it must have been
thought of before, but I have not been able to find any references.
We shall refer to it as the parallel ``self-labeling'' algorithm since in it each site
figures out which cluster its in by itself from local information.
We begin by assigning each site $i$ a unique cluster label $S_i$.
In practice this is simply chosen as the position of that site in the lattice (that is,
on an $M \times N$ lattice a site with coordinates $k,l$ has position $k + M \times l$).
We also set the elements of a count array, whose length is the number of sites in the
lattice ($M \times N$, in this example), to one: $A(S_i) = 1; i=1,2,...,M \times N$.
At each step of the algorithm, in parallel,
every site looks at each of its $d$ neighbors (here $d=2$) in the
positive directions in turn. If the neighboring site $n$ has the same spin and if it
has a different cluster label $S_n$, then both $S_i$ and $S_n$ are set to the minimum
of the two.
In addition, the count arrays are updated such that the count array for the site whose
cluster label
is overwritten is set to zero, and the count array for the other site is
incremented by one. 
In other words,
if $S_i < S_n$ then $S_n = S_i$ and $A(S_n) = 0, A(S_i) = A(S_i)+1$; otherwise,
if $S_n < S_i$ then $S_i = S_n$ and $A(S_i) = 0, A(S_n) = A(S_n)+1$.
(If $S_i = S_n$ nothing happens.)
This step is repeated until nothing changes -- which means that all connected
sites which have the same spin are labeled with the same label and thus all the
clusters are identified.
We also end up with the count array specifying how many sites are in each labeled
cluster (zero elements indicate that there is no cluster with that label).
As a final step, the count array can be scanned in order to determine the
number of clusters (at the same time it can be compacted to make it easier to deal
with later).

We note that this ``self-labeling'' algorithm works sequentially also
(as does the ``ants in the labyrinth'' algorithm).
For parallel execution we imagine one processor per lattice site.
If there are more lattice sites than processors then we employ domain decomposition
so that the processors simply cycle through the sites in their domain sequentially,
communicating across the boundaries as usual.
There is one subtlety in executing the self-labeling algorithm in parallel --
more than one processor cannot update an element of the count array at a time,
which implies that we must lock/unlock these elements.
Of course, the mechanism for this depends on which breed of parallel computer
we implement our algorithm.
On MIMD machines with no host (like the BBN Butterfly, the Sequent Balance and
Symmetry, and the Encore Multimax) we declare the count array as a shared variable
and have locks (using semaphores or fetch-and-add, for example) for each of its
elements.
On MIMD machines with a host computer (namely, hypercubes) we may choose to
store the count array on the host so that processors update its elements
by sending the host messages which are processed sequentially thus avoiding the
problem.
On SIMD machines with no host (like the AMT DAP) we must get around the problem with
extra code -- each processor notes which element of the count array it wishes to update,
then all these requests are combined (efficiently, along rows and columns in the DAP)
and the array is updated $once$.
On SIMD machines with a host (for example the TMC Connection Machine) we can
store the count array on the host and allow the processors to update it in parallel,
making use of the ``as-if-serial-rule'' to ensure that the update actually happens
as if in some sequential order and therefore gives the correct result.
Of course we can avoid this problem with the count array altogether on any
SIMD computer by getting rid of it and simply using the ``where'' and
``sum'' functions to mask off and count the clusters after they have been
built.

\bigskip
\line{\bf 7. Conclusions\hfill}

We have historically reviewed
numerical computations of discrete q-state Potts and continuous O(N) spin models,
demonstrating that the standard local update
Metropolis Monte Carlo algorithm has a dynamical critical exponent $z \approx 2$.
Next the over-relaxed algorithm was described: it allows $z \approx 1$.
Then the cluster update algorithms which can yield $z \approx 0$ were explained.
Finally we have discussed how one identifies clusters both sequentially and in
parallel.

In conclusion, the future looks bright for spin model computations -- yes,
they will all be done again! -- using cluster algorithms.

\bigskip
\line{\bf Acknowledgements\hfill}

I would like to thank Paul Coddington and Rajan Gupta for useful discussions.

\listrefs\vfill
\end
