
\section{Spin Models}
\label{spins}

\subsection{Introduction}

Spin models\index{Spin models} are simple statistical models of real systems, 
which exhibit the same behavior and, hence, provide an understanding of the 
physical mechanisms involved.  Despite their apparent simplicity, most of 
these models are not exactly soluble by present theoretical methods.  Hence, 
computer simulation is used.  Usually, one is interested in the behavior of 
the system at a phase transition; the computer simulation reveals where the 
phase boundaries are, what the phases on either side are, and how the 
properties of the system change across the phase transition.  There are two 
varieties of spins\@: discrete or continuously valued.  In both cases, the 
spin variables are put on the sites of the lattice and only interact with 
their nearest neighbors.  The partition function for a spin model is
\begin{equation}
Z = \int Ds \exp(S),
\end{equation}
with the action being of the form
\begin{equation}
S = \beta \sum_{\langle ij\rangle} ( 1 - s_i.s_j ),
\end{equation}
where $\langle ij\rangle$ denotes nearest neighbors.  A great deal of work 
has been done over the years in finding good algorithms for computer 
simulations of spin models; recently some new, much better, algorithms have 
been discovered.  These so-called cluster algorithms are described in
detail in Section~\ref{cluster}.  Here, we shall use them to perform
large-scale Monte Carlo simulations of several spin models---both
discrete and continuous.

\subsection{Ising Model}

The Ising model\index{Ising model} is the simplest model for ferromagnetism 
that predicts phase transitions and critical phenomena.  This model, 
introduced by Lenz in 1920 \cite{Lenz:20a}, was solved in one dimension by 
Ising in 1925 \cite{Ising:25a}, and in two dimensions by Onsager in 1944 
\cite{Onsager:44a}.  However, it has not been solved analytically in three 
dimensions, so Monte Carlo computer simulation methods have been one of the 
methods used to obtain numerical solutions.  One of the best available 
techniques for this is the Monte Carlo Renormalization Group\index{Monte
Carlo Renormalization Group} (MCRG) method \cite{Wilson:80a}, 
\cite{Swendsen:79a}.  The Ising model exhibits a second-order phase 
transition in $d=3$ dimensions at a critical temperature $T_c$.  As $T$ 
approaches $T_c$, the correlation length $\xi$ diverges as a power law with 
critical exponent $\nu$\@:  
\begin{equation}
\xi = \xi_0(T/T_c-1)^{-\nu}
\end{equation}
and the pair correlation function $G(r)$ at $T=T_c$ falls off to zero with
distance $r$ as a power law defining the critical exponent\index{critical 
exponents} $\eta$\@:
\begin{equation}
G(r) = G_0 r^{-(d-2+\eta)}.
\end{equation}
$T_c$, $\nu$ and $\eta$ determine the critical behavior of the 3-D Ising model
and it is their values we wish to determine using MCRG.

In 1984, this was done by Pawley, Swendsen, Wallace and Wilson
\cite{Pawley:84a} in Edinburgh on the ICL DAP computer with high statistics.
They ran on four lattice sizes---$8^3$, $16^3$, $32^3$ and $64^3$---measuring 
seven even and six odd spin operators.  We are essentially repeating their 
calculation on the new AMT DAP computer.  Why should we do this?  First, to 
investigate finite size effects---we have run on the biggest lattice used by 
Edinburgh, $64^3$, and on a bigger one, $128^3$\@.  Second, to investigate 
truncation effects---qualitatively the more operators we measure for MCRG, 
the better, so we have included 53 even and 46 odd operators.  Third, we are 
making use of the new cluster-updating algorithm due to Swendsen and Wang 
\cite{Swendsen:87a}, implemented according to Wolff \cite{Wolff:89b}.  
Fourth, we would like to try to measure another critical exponent more
accurately---the correction-to-scaling exponent $\omega$, which plays an
important role in the analysis.

In order to calculate the quantities of interest using MCRG, one must 
evaluate the spin operators $S_\alpha$.  In \cite{Pawley:84a}, the 
calculation was restricted to seven even spin operators and six odd; we are 
evaluating 53 and 46, respectively.  Specifically, we have decided to 
evaluate the most important operators in a $3\times 3\times 3$ cube 
\cite{Baillie:88h}.  To determine the critical coupling (or inverse 
temperature), $K_c$, one performs independent MC simulations on a large 
lattice $L$ of size $M^3$ and on smaller lattices $S$ of size $(M/2^m)^3$, 
$m = 1,2,\,\ldots$, and compares the operators measured on the large lattice 
blocked $m$ times more than the smaller lattices.  $K = K_c$ when they are 
the same.  Since the effective lattice sizes are the same, unknown finite 
size effects should cancel.  The critical exponents, $y_a$, are obtained 
directly from the eigenvalues, $\lambda_a$, of the stability matrix, 
$T_{\alpha\beta}$, according to $\lambda_a = 2^{y_a}$.  In particular, the 
leading eigenvalue $y_t$ of $T_{\alpha\beta}$ for the even $S_\alpha$ gives 
$\nu$ from $y_5 = 1/\nu$, and similarly $y_h = (d+2-\eta)/2$ from the odd 
eigenvalue $y_h$ of $T_{\alpha\beta}$.

The Distributed Array Processor (DAP) is a SIMD computer consisting of $D 
\times D$ bit-serial processing elements (PEs) configured as a cyclic 
two-dimensional grid with nearest-neighbor connectivity.  The Ising model 
computer simulation is well-suited to such a machine since the spins can be 
represented as single-bit (logical) variables.  In 3-D, the system of spins 
is configured as an $M \times M \times M$ simple cubic lattice, which is 
``crinkle mapped'' onto the $D \times D$ DAP by storing $N \times N$ pieces 
of each of $M$ planes in each PE\@: $M \times M \times M = M \times (N\times 
D) \times (N\times D)$, with $N = M/D$.  Our Monte Carlo simulation uses a 
hybrid algorithm in which each sweep consists of 10 standard Metropolis 
\cite{Metropolis:53a} spin-updates followed by one cluster-update using 
Wolff's single-cluster variant of the Swendsen and Wang algorithm.  On the 
$128^3$ lattice, the autocorrelation time of the magnetization reduces from 
$73 \pm 2$ sets of $100$ sweeps for Metropolis alone to $5.0 \pm 0.2$ sets of 
$10$ Metropolis plus one cluster-update for the hybrid algorithm.  In order 
to measure the spin operators, $S_\alpha$, the DAP code simply histograms the 
spin configurations so that an analysis program can later pick out each 
particular spin operator using a look-up table.  Currently, the code requires 
the same time to do one histogram measurement, one Wolff 
single-cluster-update or $100$ Metropolis updates.  Therefore, our hybrid of 
$10$ Metropolis plus one cluster-update takes about the same time as a 
measurement.  On a DAP~510, this hybrid update takes on average 127
secs (13.5 secs) for the $128^3$ ($64^3$) lattices, so that the total
time used for this calculation is roughly 11,000 hours.  We have performed
simulations on $64^3$ and $128^3$ lattices at two values of the
coupling\@:  $K_{nn} = 0.221654$ (Edinburgh's best estimate) and 
$K_{nn} = 0.221644$.  We accumulated $\rm 140\,K$ measurements for each
of the $128^3$ simulations and $\rm 100\,K$ for the $64^3$.  For error
analysis, this is divided into bins of $\rm 10\,K$ measurements.

In analyzing our results, the first thing we have to decide is the order in 
which to arrange our 53 even and 46 odd spin operators.  Naively, they can be 
arranged in order of increasing total distance between the spins 
\cite{Baillie:88h} (as was done in \cite{Pawley:84a}).  However, the ranking 
of a spin operator is determined physically by how much it contributes to 
the energy of the system.  Thus, we did our analysis initially with the 
operators in the naive order to calculate their energies, then subsequently 
we used the ``physical'' order dictated by these energies.  This physical 
order of the first 20 even operators is shown in Figure~\ref{fig1.spins} with
six of Edinburgh's operators indicated; the seventh Edinburgh operator (E-6) 
is our 21st.  This order is important in assessing the systematic effects of 
truncation, as we are going to analyze our data as a function of the 
number of operators included.  Specifically, we successively diagonalize the 
$1\times 1$, $2\times 2$, $\ldots$, (for even, $46\times 46$ for odd) 
stability matrix $T_{\alpha\beta}$ to obtain its eigenvalues and, thus, the 
critical exponents.

\begin{figure}
\vspace*{2.5in}
\caption{Our Order for Even Spin Operators}
\label{fig1.spins}
\end{figure}

We present our results in terms of the eigenvalues of the even and odd parts 
of $T_{\alpha\beta}$.  The leading even eigenvalue on the first four blocking
levels starting from the $64^3$ lattice is plotted against the number of 
operators included in the analysis in Figure~\ref{fig2.spins}, and on the 
first five blocking levels starting from the $128^3$ lattice in 
Figure~\ref{fig3.spins}.  Similarly, the leading odd eigenvalues for $64^3$ 
and $128^3$ lattices are shown in Figure~\ref{fig4.spins} and
Figure~\ref{fig5.spins} respectively.  First of all, note that there are 
significant truncation effects---the value of the eigenvalues do not settle 
down until at least thirty and perhaps forty operators are included.  We note 
also that our value agrees with Edinburgh's when around seven operators are
included---this is a significant verification that the two calculations are
consistent.  With most or all of the operators included, our values on the 
two different lattice sizes agree, and the agreement improves with increasing 
blocking level.  Thus, we feel that we have overcome the finite size effects 
so that a $64^3$ lattice is just large enough.  However, the advantage in 
going to $128^3$ is obvious in Figure~\ref{fig3.spins} and
Figure~\ref{fig5.spins}\@:  there, we can perform one more blocking, which 
reveals that the results on the fourth and fifth blocking levels are 
consistent.  This means that we have eliminated most of the transient effects 
near the fixed point in the MCRG procedure.  We also see that the main
limitation of our calculation is statistics--the error bars are still
rather large for the highest blocking level.

Now in order to obtain values for $\nu$ and $\eta$, we must extrapolate
our results from a finite number of blocking levels to an infinite
number.  This is done by fitting the corresponding eigenvalues
$\lambda_t$ and $\lambda_h$ according to
\begin{equation}
\lambda(n) = \lambda^{\ast} + a2^{-\omega n}
\end{equation}
where $\lambda^{\ast}$ is the extrapolated value and $\omega$ is the
correction-to-scaling exponent.  Therefore, we first need to calculate
$\omega$, which comes directly from the second leading even
eigenvalue\@: $\omega = -y_{t,2}$.  Our best estimate is in the interval
$\omega = 0.80$--0.85, and we use the value 0.85 for the purpose of
extrapolation, since it gives the best fits.  The final results are
$\nu = 0.624\pm 0.001\pm 0.002$, $\eta = 0.0262\pm 0.0002\pm 0.0030$,
where the first errors are statistical and the second errors are
estimates of the systematic error coming from the uncertainty in
$\omega$.

\begin{figure}
\vspace*{2.0in}
\caption{Leading Even Eigenvalue on $64^3$ Lattice}
\label{fig2.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Even Eigenvalue on $128^3$ Lattice}
\label{fig3.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Odd Eigenvalue on $64^3$ Lattice}
\label{fig4.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Odd Eigenvalue on $128^3$ Lattice}
\label{fig5.spins}
\end{figure}

Finally, perhaps the most important number, because it can be determined the 
most accurately, is $K_{nn}^c$.  By comparing the fifth blocking level on the 
$128^3$ lattice to the fourth on the $64^3$ lattice for both coupling
values and taking a weighted mean, we obtain $K_{nn}^c = 0.221652\pm
0.000003\pm 0.000001$, where again the first error is statistical and
the second is systematic.

Thus, MCRG calculations give us very accurate values for the three critical 
parameters $K_{nn}^c$, $\eta$, and $\nu$, and give a reasonable estimate
for $\omega$.  Each parameter is obtained independently and directly
from the data.  We have shown that truncation and finite size errors at
all but the highest blocking level have been reduced to below the
statistical errors.  Future high statistics simulations on $256^3$
lattices will significantly reduce the remaining errors and allow us to
determine the exponents very accurately.

\subsection{Potts Model}
\label{potts}

The $q$-state Potts model\index{Potts model} \cite{Potts:52a} consists of a 
lattice of spins $\sigma_{i}$, which can take $q$ different values, and whose 
Hamiltonian is
\begin{equation}
H = K \sum_{\langle i,\,j\rangle} \delta_{\sigma_{i}\sigma_{j}} .
\end{equation}
For $q=2$, this is equivalent to the Ising model.  The Potts model is thus a 
simple extension of the Ising model; however, it has a much richer phase 
structure, which makes it an important testing ground for new theories and
algorithms in the study of critical phenomena \cite{Wu:82a}. 

Monte Carlo simulations of Potts models have traditionally used local 
algorithms such as that of Metropolis, {\em et al.} \cite{Metropolis:53a},
however, these algorithms have the major drawback that near a phase 
transition the autocorrelation time (the number of sweeps needed to generate 
a statistically independent configuration) increases approximately as $L^2$, 
where $L$ is the linear size of the lattice.  New algorithms have recently 
been developed that dramatically reduce this ``critical slowing
down''\index{critical slowing down} by 
updating clusters of spins at a time (these algorithms are described in 
Section~\ref{cluster}).  The original cluster algorithm of Swendsen and Wang 
was implemented for the Potts model \cite{Swendsen:87a}, and there is a lot 
of interest in how well cluster algorithms perform for this model.  At 
present, there are very few theoretical results known about cluster 
algorithms, and theoretical advances are most likely to come from first 
studying the simplest possible models.

We have made a high statistics study of the SW algorithm and the single 
cluster Wolff algorithm \cite{Wolff:89b}, as well as a number of variants of 
these algorithms, for the $q=2$ and $q=3$ Potts models in two dimensions 
\cite{Baillie:89m}.  We measured the autocorrelation time $\tau$ in the 
energy (a local operator) and the magnetization (a global one) on lattice 
sizes from $8^2$ to $512^2$.   About 10 million sweeps were required for each 
lattice size in order to measure autocorrelation times to within about 1\%.
From these values, we can extract the dynamic critical exponent $z$, given 
by $\tau(L) \sim L^{z}$, where $\tau(L)$ is measured at the infinite volume 
critical point (which is known exactly for the $2D$ Potts model). 

The simulations were performed on a number of different parallel computers. 
For lattice sizes of $128^2$ or less, it is possible to run independent
simulations on each processor of a parallel machine, enabling us to obtain
100\% efficiency by running 10 or 20 runs for each lattice size in parallel.
These calculations were done using a 32-processor Meiko Computing Surface,
a 20-processor Sequent Symmetry, a 20-processor Encore Multimax, and
a 96-processor BBN GP1000 Butterfly, as well as a network of SUN workstations,
and took approximately 15,000 processor-hours.  For the largest lattice sizes, 
$256^2$ and $512^2$, a parallel cluster algorithm was required, due to the 
large amount of calculation (and memory) required.  We have used the 
self-labeling algorithm described in Section~\ref{cluster}, which gives 
fairly good efficiencies of about 70\% on the machines we have used (an 
nCUBE-1 and a Symult S2010), by doing multiple runs of 32 nodes each for the 
$256^2$ lattice, and 64 nodes for $512^2$.  Since this problem does not
vectorize, using all 512 nodes of the nCUBE gives a performance approximately 
five times that of a single processor CRAY~X-MP, while all 192 nodes of the 
Symult is equivalent to about six CRAYs.  The calculations on these machines 
have so far taken about 1000 hours.

Results for the autocorrelation times of the energy for the Wolff and SW 
algorithms are shown in Figure~\ref{fig6.spins} for $q=2$ and
Figure~\ref{fig7.spins} for $q=3$.  As can be seen, the Wolff algorithm has 
smaller autocorrelation times than SW.  However, the dynamical critical 
exponents for the two algorithms appear to be identical, being approximately 
$0.25(1)$ and $0.57(1)$ for $q=2$ and ${q=3}$ respectively (shown as
straight lines in Figure~\ref{fig6.spins} and Figure~\ref{fig7.spins}), 
compared to values of approximately two for the standard Metropolis 
algorithm. 

\begin{figure}
\vspace*{2.0in}
\caption{Energy Autocorrelation Times, $q=2$}
\label{fig6.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Energy Autocorrelation Times, $q=3$}
\label{fig7.spins}
\end{figure}

Burkitt and Heermann \cite{Heermann:90a} have suggested that the increase in 
the autocorrelation time is a logarithmic one, rather than a power law for 
the $q=2$ case (the Ising model), {\em i.e.,\/} $z = 0$.  Fits to this
are shown as dotted lines in Figure~\ref{fig6.spins}.  These have
smaller $\chi^2$ values than the power law fits, favoring logarithmic
behavior.  However, it is very difficult to distinguish between a logarithm 
and a small power even on lattices as large as $512^2$.  In any case, the 
performance of the cluster algorithms for the Potts model is quite 
extraordinary, with autocorrelation times for the $256^2$ lattice, which are 
hundreds of times smaller than for the Metropolis algorithm.  In the future, 
we hope to use the cluster algorithms to perform greatly improved Monte Carlo 
simulations of various Potts models, to study their critical behavior.

\subsection{XY Model}

The XY (or O(2)) model\index{XY model} consists of a set of continuous valued 
spins regularly arranged on a two-dimensional square lattice.  15 years ago, 
Kosterlitz and Thouless\index{Kosterlitz and Thouless} (KT) predicted that 
this system would undergo a phase transition\index{transitions} as one 
changed from a low temperature spin wave phase to a high 
temperature phase with unbound vortices.  KT  predicted an approximate 
transition temperature, $T_c$, and the following unusual exponential 
singularity in the correlation length and magnetic susceptibility\@:
\begin{equation}
\xi = a_{\xi} e^{b_{\xi}t^{-\nu}},\quad\quad \chi = 
a_\chi e^{b_\chi t^{-\nu}},
\end{equation}
with
\begin{equation}
\nu = \frac{1}{2},\quad\quad \eta = \frac{1}{4},
\end{equation}
where $t = T-T_c$ and the correlation function exponent $\eta$ is defined by 
the relation $\chi = c\, \xi^{2-\eta}$.  

Our simulation \cite{Gupta:88a} was done on the 128-node FPS (Floating Point 
Systems) T-Series hypercube at Los Alamos.  FPS software allowed the use of C 
with a similar software model (communication implemented by subroutine call) 
to that used on the hypercubes at Caltech.  Each FPS node is built around 
Weitek floating point units, and we achieved $\rm 2\,Mflops$ per node in this 
application.  The total machine ran at $\rm 1/4\,Gflops$, or at about twice 
the performance of one processor of a CRAY~X-MP for this application.  We use 
a $1-D$ torus topology for communications, with each node processing a 
fraction of the rows.  Each row is divided into red/black alternating sites of
spins and the vector loop is over a given color.  This gives a natural data 
structure of ($n \times 128$) words for lattices of size $(n \times 256)$.  
The internode communications, in both lattice update and in measurement of 
observables, can be done asynchronously and are a negligible overhead.

Previous numerical work was unable to confirm the KT theory, due to limited 
statistics and small lattices.  Our high-statistics simulations are done on 
$64^2$, $128^2$, $256^2$, and $512^2$ lattices using a combination of 
over-relaxed and Metropolis algorithms which decorrelates as $\tau \approx 
0.15 \xi^{1.2}$.  (For comparison, a Metropolis algorithm decorrelates as 
$\tau \approx 5 \xi^2$.)  Each configuration represents $N_{or}$ over-relaxed 
sweeps through the lattice followed by $N_{met}$ Metropolis sweeps.  
Measurement of observables is made on every configuration.  The over-relaxed 
algorithm consists of reflecting the spin at a given site about $\Sigma$ 
where $\Sigma$ is the sum of the nearest-neighbor spins, {\em i.e.},
\begin{equation}
s_{new} \ \ = \ \ \Sigma \ s_{old}^\dagger \ \Sigma.
\end{equation}
This implementation \cite{Creutz:87a} \cite{Brown:87a} of the over-relaxed 
algorithm is microcanonical, and it reduces critical slowing down even though 
it is a local algorithm.  The ``hit'' elements for the Metropolis algorithm
are generated as $e^{i\theta}$, where $\theta$ is a uniform random number
in the interval $\lbrack - \delta, \delta\rbrack$, and $\delta$ is adjusted to 
give an acceptance rate of 50--60\%.  The Metropolis hits make the algorithm 
ergodic, but their effectiveness is limited to local changes in the energy.
In Figure~\ref{fig8.spins}, we show the autocorrelation time $\tau$ vs. the 
correlation length $\xi$\@; for $N_{or} = 8$,  $N_{met} = 2$ we extract $\tau 
= 0.15\xi^{1.48}$ and for $N_{or} = 15$, $N_{met} = 2$ we get $\tau = 
0.15\xi^{1.20}$.

\begin{figure}
\vspace*{2.0in}
\caption{Autocorrelation Times for the XY Model}
\label{fig8.spins}
\end{figure}

We ran at 14 temperatures near the phase transition and made unconstrained 
fits to all 14 data points (four parameter fits according to equations~(XY)),
for both the correlation length (Figure~\ref{fig9.spins}) and susceptibility 
(Figure~\ref{fig10.spins}).  The key to the interpretation of the data is the 
fits.  We find that fitting programs (e.g., MINUIT, SLAC) move incredibly
slowly towards the true minimum from certain points (which we label spurious
minima), which, unfortunately, are the attractors for most starting points.
We found three such spurious minima (KT1-3) and the true minimum KT4, as 
listed in Table~\ref{tab1.spins}.

\begin{figure}
\vspace*{2.0in}
\caption{Correlation Length for the XY Model}
\label{fig9.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Susceptibility for the XY Model}
\label{fig10.spins}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|l|}
\hline
&\multicolumn{1}{|c|}{KT1}&\multicolumn{1}{c|}{KT2}&\multicolumn{1}{c|}{KT3}&
\multicolumn{1}{c|}{KT4}
\\
\hline
\medStrut
$a_{\chi}$      & "2.155 & "0.270 & 0.784 & 0.2727(3)  \\
$b_{\chi}$      & "1.088 & "2.809 & 1.844 & 2.7734(5)  \\
$T_c$           & "0.701 & "0.911 & 0.867 & 0.8987(1)  \\
$\nu$           & "1.460 & "0.472 & 0.701 & 0.4995(11) \\
$\chi^2/_{dof}$ & 73.0"" & 13.2"" & 1.82" & 1.70       \\
$a_{\xi}$       & "0.708 & "0.228 & 0.441 & 0.2280(3)  \\
$b_{\xi}$       & "0.740 & "1.693 & 1.100 & 1.6725(8)  \\
$T_c$           & "0.689 & "0.910 & 0.867 & 0.8967(2)  \\
$\nu$           & "1.412 & "0.471 & 0.700 & 0.5007(3)  \\
$\chi^2/_{dof}$ & "7.94" & "2.62" & 0.89" & 0.37       \\
$\eta$          & "0.53" & "0.34" & 0.324 & 0.34       \\
\hline
\end{tabular}
\caption{Results of the XY model fits\@: (a) $\chi$ in T, and  (b) $\xi$ 
in T assuming the KT form.  The fits KT1-3 are pseudo-minima while KT4 is 
the true minimum.  All data points are included in the fits and we give the 
$\chi^2/_{dof}$ for each fit and an estimate of the exponent $\eta$.}
\label{tab1.spins}
\end{center}
\end{table}

Thus, our data was found to be in excellent agreement with the KT theory and, 
in fact, this study provides the first direct measurement of $\nu$ from both 
$\chi$ and $\xi$ data that is consistent with the KT predictions.

\subsection{O(3) Model}\index{O(3) Model}

The XY model is the simplest O(N) model, having $N=2$; the O(N) model being a 
set of rotors (N-component continuous valued spins) on an N-sphere.  For $N 
\ge 3$, this model is asymptotically free \cite{Polyakov:75a}, and for $N=3$, 
there exist instanton solutions.  Some of these properties are analogous to 
those of gauge theories in four dimensions; hence, these models are 
interesting.  In particular, the O(3) model\index{O(3) model} in two 
dimensions should shed 
some light on the asymptotic freedom of QCD (SU(3)) in four dimensions.
The predictions of the renormalization group for the susceptibility $\chi$
and inverse correlation length ({\em i.e.,\/} mass gap) $m$ in the O(3) model 
are \cite{Brezin:76a}
\begin{equation}
\chi = C \beta^{-4} \exp\,(4\pi\beta )\, (1 + O(1/\beta))
\end{equation}
and
\begin{equation}
m = \frac{1}{\xi} = C^\prime \beta \exp\,(-2\pi\beta )\,(1 + O(1/\beta)) ,
\end{equation}
respectively.
If $m$ and $\chi$ vary according to these equations, without the correction of
order $1/\beta$, they are said to follow asymptotic scaling.  Previous work 
was able to confirm that this picture is qualitatively correct, but was not 
able to probe deep enough in the area of large correlation lengths to obtain 
good agreement. 
 
The combination of the over-relaxed algorithm and the computational power of 
the FPS T-Series allowed us to simulate lattices of size up to $1024^2$.  We
were thus, able to simulate at coupling constants that correspond to 
correlation lengths up to $300$, on lattices where finite size effects are 
negligible.  We were also able to gather large statistics and, thus, obtain 
small statistical errors.  Our simulation is in good agreement with the 
similar cluster calculations \cite{Wolff:89b} \cite{Wolff:90a}.  Thus, we 
have validated and extended these results in a regime where our algorithm is 
the only known alternative to clustering.

We have made extensive runs at 10 values of the coupling constant.  At the 
lowest $\beta$, several hundred thousand sweeps were collected while for the 
largest values of $\beta$, between $50,000$ and $100,000$ sweeps were made.  
Each sweep consists of between 10 iterations through the lattice at the 
former end and 150 iterations at the latter.  The statistics we have gathered 
are equivalent to about 200 days, using the full 128-node FPS machine.

Our results for the correlation length and susceptibility for each coupling
constant $\beta$ and lattice size, are shown in Table~\ref{tab2.spins}.  The 
autocorrelation times are also shown.  The quantities measured on different 
sized lattices at the same $\beta$ agree, showing  that the infinite volume 
limit has been reached.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{Beta}&\multicolumn{1}{c|}{L,T}&
\multicolumn{1}{c|}{$\tau$}&\multicolumn{1}{c|}{$N_{or}$}&
\multicolumn{1}{c|}{$\chi$}&\multicolumn{1}{c|}{$\xi$}
\\
\hline
\medStrut
1.50 & "256     & 1.3 & "12 & $""176.4 \pm ""0.2$ & "11.1 \\
1.60 & "256     & 2.5 & "12 & $""448.4 \pm ""0.7$ & "19.0 \\
1.70 & "256     & 7.0 & "12 & $"1263\phantom{.}"\pm ""3$\phantom{.}" & "34.4 \\
     & "512     & 1.8 & "35 & $"1264\phantom{.}"\pm ""3$\phantom{.}" & "34.4 \\
1.75 & "768     & 1.8 & "50 & $"2197\phantom{.}"\pm "15$\phantom{.}" & "47.2 \\
1.80 & "512     & 3.5 & "40 & $"3850\phantom{.}"\pm "10$\phantom{.}" & "64.7 \\
     & "768     & 3.0 & "45 & $"3820\phantom{.}"\pm "20$\phantom{.}" & "64.5 \\
1.85 & "768     & 2.3 & "80 & $"6730\phantom{.}"\pm "25$\phantom{.}" & "88.7 \\
1.90 & 1024, 512& 2.0 & 150 & $11600\phantom{.}"\pm "60$\phantom{.}" & 121.5 \\
     & 1024     & 1.5 & 200 & $11870\phantom{.}"\pm "60$\phantom{.}" & 122.7 \\
1.95 & 1024     & 1.7 & 200 & $20640\phantom{.}"\pm 310$\phantom{.}" & 164.8 \\
2.00 & 1024     & 1.7 & 250 & $35100\phantom{.}"\pm 400$\phantom{.}" & 224.3 \\
2.05 & 1024     & 1.7 & 300 & $56220\phantom{.}"\pm 550$\phantom{.}" & 295.6 \\
\hline
\end{tabular}
\caption{Coupling constant, lattice size, autocorrelation time, number of 
over-relaxed sweeps, susceptibility and correlation length for the O(3) 
model.}
\label{tab2.spins}
\end{center}
\end{table}

To compare the behavior of the correlation length and susceptibility
with the asymptotic scaling predictions, we use the ``correlation length
defect'' $\delta_{\xi}$ and ``susceptibility defect'' $\delta_{\chi}$
which are defined as follows\@: $\delta_{\xi} = \beta (e^{-2\pi\beta})\xi$, 
$\delta_{\chi} = \beta^4 (e^{-4\pi\beta})\chi$ so that asymptotic
scaling is seen if $\delta_{\xi}$, $\delta_{\chi}$ go to constants as
$\beta\rightarrow\infty$.  These defects are shown in
Figure~\ref{fig11.spins} and Figure~\ref{fig12.spins} respectively.  It
is clear that asymptotic scaling does not set in for $\beta < 2$, but it
is not possible to drawn a clear conclusions for $\beta > 2$---though
the trends of the last two or three points may be toward constant
behvaior.

\begin{figure}
\vspace*{2.0in}
\caption{Correlation Length Defect vs. the Coupling Constant for the O(3) 
Model}
\label{fig11.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Susceptibility Defect vs. the Coupling Constant for the O(3) Model}
\label{fig12.spins}
\end{figure}

We gauged the speed of the algorithm in producing statistically independent 
configurations by measuring the autocorrelation time $\tau$.  We used this to 
estimate the dynamical critical exponent $z$, which is defined by $\tau = c
\cdot\xi^z$.  For constant $N_{or}=12$, our fits give $z=1.33(1)$.  However, 
we discovered that by increasing $N_{or}$ in rough proportion to $\xi$, we 
can improve the performance of the algorithm significantly.  To compare the
speed of decorrelation between runs with different $N_{or}$, we define a new 
quantity, which we call `effort', $e =  N_{or} \times\tau$.  This measures 
the computational effort expended to obtain a decorrelated configuration.  We 
define a new exponent $z^\prime$  from $e\sim\xi^{z^\prime}$, where $N_{or}$ 
is chosen to keep $\tau$ constant.  We also found that the behavior of
the decorrelation time can be approximated over a good range by
\begin{equation}
\tau  = C^{\prime\prime}\cdot\xi^z \cdot N_{or}^{-z/z^\prime}.
\end{equation}
A fit to the set of points ($N_{or}$, $\xi$, $\tau > 1.0$) gives
$z = 1.301\pm 0.012$, $z^{\prime} = 1.079\pm 0.010$.  Thus, $z^{\prime}$ is
significantly lower than $z$.  Figure~\ref{fig13.spins} shows $\tau$ vs
$N_{or}$ with the fits shown as solid lines.

\begin{figure}
\vspace*{2.0in}
\caption{Decorrelation Time $\tau$ vs. Number of Overrelations Sweeps
$N_{or}$ for Different Values of $\beta$}
\label{fig13.spins}
\end{figure}

