\documentstyle[aps,epsf,preprint]{revtex}
\def\boldmu{\mbox{\boldmath $\mu$}}
\def\hatmu{\hat \mu}
\def\H0{{\bf H}_0}
\def\boldH{{\bf H}}
\def\boldr{{\bf r}}

\begin{document}
\draft

\title{Cluster-based Monte Carlo Simulation of Ferrofluids}

\author{S.W. Davis, W. McCausland, H.C. McGahagan, C.T. Tanaka and M. Widom}
\address{Department of
Physics, Carnegie-Mellon University\\ Pittsburgh, Pa. 15213}

\date{\today}
\maketitle

\begin{abstract}
We demonstrate a Monte Carlo algorithm for efficiently simulating
ferrofluids. By identifying particle clusters and evolving them as
single units, we reduce correlation times by more than 2 orders of
magnitude. This method enables accurate calculations of ferrofluid
thermodynamics in the limit of strong magnetic coupling that would be
impossible by conventional means. We apply the method to study
magnetic anisotropy in dilute thin films.
\end{abstract}

\pacs{PACS: 75.50.Mm, 61.20.Ja, 75.70.Ak, 82.70.Dd}

\section{Introduction}
Ferrofluids~\cite{Rosensweig} are colloidal suspensions of
nanometer-scale magnetic particles. When the magnetic interactions
between particles are stronger than thermal energy, the particles
aggregate into chains~\cite{deGennes}, with the magnetic moment of
each particle pointing towards a neighbor along the chain (see
figure~\ref{conf}).  Theoretical analysis of liquid state
thermodynamics in this interesting state is difficult. Weak magnetic
coupling expansions are not accurate because particle pairs are
tightly bound~\cite{smalllambda}. Virial
expansions~\cite{virial,WZ,thesis} in low particle density are not
adequate either, because chain formation is inherently a many-particle
correlation.  While there is hope for analytic progress using concepts
of ``living polymers''~\cite{WZ,living}, it is important to test such
approximate theories with computer simulation.

Traditional Monte Computer simulation, likewise, is difficult when
particle chains form~\cite{chainsim}.  The problem is that the strong
attraction of particles within a chain require very small particle
steps to achieve a reasonable ratio of accepted steps. Large-scale
translations, rotations and undulations of the chains, which are an
important factor in determining such important thermodynamic
quantities as the magnetic susceptibility, are poorly simulated using
only independent small displacements of each particle. Evidence of the
difficulty is evident in table~\ref{std_data}, showing the growth of
correlation times as temperature drops.

This paper explores cluster-based enhancements to the traditional
Monte Carlo technique that identify particle clusters and evolve them
as entire units. Similar ideas have previously been employed in
simulations of strongly interacting ionic fluids~\cite{Ork}.
We describe how to combine traditional
single-particle steps with multi-particle steps, while maintaining
detailed balance. We employ cluster moves for speeding, to speed
magnetization reversal and chain translation. Our results show that
cluster-based algorithms are of substantial utility, speeding
equilibration of internal energy by a factor of 2 (in CPU time), and
speeding equilibration of magnetization by a factor of 360. As a
result we achieve quick, accurate calculations of susceptibility that
are difficult with conventional techniques.

After establishing the utility of cluster-based Monte Carlo
simulations, we apply the method to examine magnetic anisotropy in
ferrofluid thin films. Macroscopic thin films exhibit effective
anisotropy due to demagnetization effects. Microscopic thin films
exhibit dramatically enhanced anisotropy (see table~\ref{std_data})
due to confinement of clusters within a plane.

\section{Model and Algorithms}

We model the ferrofluid as a gas of dipolar hard spheres~\cite{DHS}.
Each sphere has diameter $a$ and magnetic moment $\mu$. One defines a
reduced temperature $T^{\star} = k_B T a^3/\mu^2$ and a reduced
density $\rho^{\star} = \rho a^3$. The dipole moment of particle $i$
is a three dimensional vector $\boldmu_i$, and the dipole interaction
falls off as $1/r^3$. To avoid the difficulty of summing the dipole
interaction in three-dimensional space~\cite{BC}, we choose to confine
our particle positions $\boldr_i$ to a two-dimensional plane.  Ewald
summation and special boundary treatment is therefore not needed. We
employ the ``minimum image'' convention for periodic boundary
conditions so that each particle interacts with the nearest image of
every other particle. In addition to simplifying our treatment of the
long-range dipolar interaction, this pseudo-2d geometry allows us to
test a theory of magnetic anisotropy in thin films~\cite{WZ}. A
related model, in which the dipole moment also is constrained to the
plain, has been examined by Weis~\cite{Weis}.

We use the Metropolis Monte Carlo method~\cite{Metropolis}, following
a Markov Chain through configuration space. Each Monte Carlo step is
accepted with probability 1 if the resulting energy change $\Delta E <
0$, and accepted with probability $\exp{(-\Delta E/k_B T)}$ otherwise.
This method produces the correct equilibrium ensemble in the limit of
infinite simulation time when certain conditions are met: the
simulation must be ergodic, so all configurations can be reached in
principle; attempted steps obey detailed balance, so no bias other
than energy influences the distribution of configurations.

These conditions are met by a conventional algorithm consisting of two
types of Monte Carlo step: randomly choose a particle, and attempt to
move it through a randomly chosen small displacement; randomly choose
a particle and attempt to rotate it by a randomly chosen angle about a
randomly chosen axis. The difficulty with this technique is the rate
at which thermodynamic properties converge to their long time limit.
Consider a particle within a chain at a temperature such that
$T^{\star} << 1$. In order for steps and rotations to be accepted with
probability near 0.5 (as is the case for the runs in
table~\ref{std_data}) their size must be very small at low
temperatures: step sizes fall off proportionally to $T^{\star}$, while
rotation angles fall off like $\sqrt{T^{\star}}$. Such small
displacements and rotations do not significantly alter the positions
and orientations of chains, or the distribution of chain lengths.

In the following we shall refer to correlation times in units of Monte
Carlo steps (MCS). One MCS in the conventional algorithm denotes $N$
attempted particle moves and $N$ attempted particle rotations. i.e.
the typical time in which each particle has a chance to move and
rotate once.  Note that it takes a time of order $N$ to evaluate potential
energy for a single particle, so the CPU time for one MCS grows like
$N^2$. Algorithms exist~\cite{multipole} that reduce this exponent but
they are not efficient for the numbers of particles considered in this
paper. For $N=100$ particles one conventional MCS takes time $t
\approx 0.2$~sec on a 50 MHz SPARC 10 workstation. We find correlation
times (see table~\ref{std_data}) in excess of $10^5$~MCS at low
temperatures, meaning it requires over five hours to generate each
independent configuration.

We supplement this conventional method with additional Monte Carlo
steps designed to accelerate motion through the configuration space.
Because we only seek thermodynamic information from our simulation, we
may exploit non-physical dynamics provided the proper equilibrium
ensemble is maintained. This possibility is an advantage that Monte
Carlo simulation may enjoy that is not possible within molecular
dynamics~\cite{MD}, another important simulation method. Thus, we
introduce a move that reverses the magnetic moment of each particle
within a chain, and a move that translates an entire cluster in a
single step.  Because interactions of a chain with the remainder of
the system are weak compared to interactions within the chain, each of
these moves will be more easily accepted than if we attempted the same
move one particle at a time.

The requirement of detailed balance dictates that the probability to
make a move, relative to the probability to make its inverse, should
depend only on the relative probability of finding the system in the
initial and final configurations. In equilibrium this equals
$\exp{(-\Delta E/k_B T)}$ where $\Delta E$ is the energy difference
between the two configurations. We choose to identify our clusters
based on proximity of particles to each other. Particles $i$ and $j$
are in proximity when $| \boldr_i - \boldr_j | < a+\epsilon$.
Typically we choose $\epsilon = 0.2a$, corresponding to a broad
minimum in the pair distribution function.

Having identified a cluster, we must ensure that any move the cluster
makes will not prevent finding the inverse move at a later time.
Particle connectivity must not be altered.  There are several ways to
enforce this. One possibility is to leave positions unchanged when
acting on a cluster. Fortunately, reversal of cluster magnetization is
an important process that needs to be accelerated but does not alter
particle connectivity. This is the process by which we accelerate
equilibration of the magnetization. A second possibility is to move
clusters subject to the condition that the new cluster position is at
least distance $\epsilon$ from any other cluster or particle. This
type of cluster move is important for the equilibration of particle
positions and it accelerates the equilibration of energy. We have
checked that our cluster Monte Carlo method correctly reproduces the
thermodynamic properties calculated by the conventional Monte Carlo
method to within our calculated uncertainties.

We identify clusters by comparing all particle positions pair-wise.
When a connected pair is found, cluster assignments of particles are
updated, and then the pairwise comparisons are continued. The time for
this algorithm is dominated by the $N(N-1)/2$ required comparisons
because we have no a-priori knowledge of particle proximity. More
efficient algorithms are possible at the cost of more complex
programming.  This algorithm suffices because it remains within the
bound of $N^2$ time required by the conventional method.

Having paid a cost to identify clusters, we attempt many cluster moves
and reversals to take maximum advantage of our work. Thus we attempt
$N$ translations and $N$ reversals of randomly chosen clusters in each
Monte Carlo step of our modified program. This is in addition to the
$N$ attempted moves and rotations of individual particles, which we
still must preform since our cluster method is not by itself ergodic.
Recall that the cluster moves are purposely non-ergodic because we
forbid cluster attachment in order to preserve detailed balance. The
time for a full Monte Carlo step still varies as $N^2$, but with a
larger coefficient.

Initial configurations for our simulations are created by placing
particles on a hexagonal lattice with dipole moments randomly
oriented. Prior to accumulating thermodynamic and correlation function
data the system is equilibrated for a time large compared to both the
duration of the initial transient in internal energy and the largest
1/e time of fluctuations in equilibrium. The system may be regarded as
well equilibrated by this time.

Figure~\ref{gEofs} displays time-correlation functions for energy
defined as
\begin{equation}
\label{eqgofs}
g_E(t) = {{1}\over{T_{sim}-t}}
		\sum_{s=1}^{T_{sim}-t} (E(s+t)-\bar{E})(E(s)-\bar{E})
\end{equation}
where $E(t)$ is the total energy of the configuration after $t$ Monte Carlo
steps and $\bar{E}$ is the mean over the entire simulation time
$T_{sim}$. Data are shown for the conventional simulation (solid
line), a run including magnetization reversals but not cluster
translations (dashed line), and a run including all types of Monte Carlo move
(dotted line). Figure~\ref{gMofs} shows the same for magnetization
$M(t)$. Both figures correspond to the density and temperature
($\rho^{\star}=0.3$, $T^{\star}=0.2$) for which figure~\ref{conf}
shows a typical configuration.

Correlation times and run times are tabulated in Table~\ref{sim_type}.
We list the $1/e$ time, rather than the rate of exponential decay
at long times, because this can be calculated more accurately. Inspection
of figures~\ref{gEofs} and~\ref{gMofs} confirms that trends in $1/e$ times
reproduce trends in conventional correlation time obtained from slopes
of the logarithm of the correlation function.

Examine first the energy data. The conventional simulation requires
2500 MCS to reach its first $1/e$ decay of $g_E$. The magnetization
reversal improves matters slightly in units of MCS, but when the
relative run time is factored in, it does little good. On the other
hand, the cluster translations greatly speed equilibration. The $1/e$
time drops by a factor 14 in units of MCS and 2.4 in CPU seconds.
Cluster moves also speed up the approach to equilibrium of
configurations initially far removed from the equilibrium state. At
densities such that the typical inter-chain distance is large compared
to the size of a chain the speedup of energy equilibration can exceed
an order of magnitude measured in CPU seconds.

The improvement is even more dramatic for the magnetization. In this
case, the magnetization reversal and cluster translation together
lower the $1/e$ time by a factor 2000 in MCS and 360 in CPU seconds.
This is because magnetization of a cluster may be entirely reversed in
a single Monte Carlo step, while it may be nearly unachievable in the
conventional approach. Acceptance ratios for magnetization reversal
not very low, about 0.15, indicating that the typical energy cost of a
cluster reversal is significantly less than the energy cost of
reversing a single particle magnetization. This is expected because
different clusters interact only weakly with each other. As a
consequence of the reduction in $\tau_M$, accurate calculation of
magnetic susceptibility becomes practical, even in the presence of
chain formation when conventional methods require absurdly long runs.

\section{magnetic Anisotropy}

With reliable values for the susceptibility in hand, what can we do
with the results? Table~\ref{std_data} reveals growing magnetic
anisotropy as temperature drops. Particle clustering enhances the
in-plane susceptibility $\chi_{xy}$ and diminishes the out-of-plane
susceptibility $\chi_z$. This behavior is also evident in
figure~\ref{fig:chi} which displays the variation of magnetic
susceptibility with density $\rho^{\star}$ at fixed temperature
$T^{\star}=0.2$.  We can explain this effect as a result of chain
formation. The following discussion presents a semi-quantitative
analysis.

Isolated, noninteracting, particles exhibit Langevin behavior with each
particle contributing 
\begin{equation}
\label{chi1}
\chi_1 = {{\mu^2}\over{3 k_B T}}
\end{equation} 
to the susceptibility in any direction. One factor of $\mu$ arises
from the coupling of each moment to an applied field, while the other
factor represents the magnetization arising from the aligned moment.
The factor of $3$ comes because each moment rotates in a 3-dimensional
space. For particles in a straight chain of length $L$, the chain
behaves like an effective particle of dipole moment $L \mu$ and
contributes $(L \mu)^2/(2 k_B T)$ to $\chi_{xy}$. The factor of 2
replaces the factor of 3 because the chain is constrained to line in
the $xy$ plane. Collecting $N$ particles into $N/L$ chains of length
$L$ the resulting susceptibility per particle is
\begin{equation}
\label{chixy}
\chi_{xy} \approx L {{\mu^2}\over{2 k_B T}}
\end{equation}
revealing the enhancement of $\chi_{xy}$ due to particle chaining.  In
contrast, the out-of-plane susceptibility $\chi_z$ is reduced by
particle clustering because particles in clusters point towards each
other in the $xy$ plane and cannot easily rotate out of this plane.
We may estimate the susceptibility per particle as
\begin{equation}
\label{chiz}
\chi_z \approx {{N_1}\over{N}} \chi_1
\end{equation}
where $N_1$ is the number of isolated particles not belonging to
clusters.

Both approximations~(\ref{chixy}) and~(\ref{chiz}) provide excellent
approximations to the data in figure~\ref{fig:chi} at low density,
below $\rho^{\star} = 0.01$. At higher densities, chain lengths begin
to grow large compared with the persistence length of the chains.
Loops and branched networks begin to form reducing the susceptibility
below the values suggested by~(\ref{chixy}). Likewise,
equation~(\ref{chiz}) breaks down when so few particles remain unbound
that the residual $\chi_z$ contributed by bound particles dominates,
and $\chi_z$ becomes independent of density. Above $\rho^{\star}
\approx 0.05$ the in-plane susceptibility $\chi_{xy}$ reaches a
maximum and begins to decrease with increasing density. This may be
related to formation of gel-like connected networks or to phase
coexistence.  Further study is needed to clarify this behavior.

The solid and dashed curves are predictions of the magnetic virial
expansion~\cite{virial,WZ,thesis}. The noninteracting Langevin value
$\chi_1$, given by equation~(\ref{chi1}), is indicated on the vertical
axis. The solid lines include the contributions of the second virial
coefficients, and the dashed lines include the contributions of the
third virial coefficients~\cite{virial,WZ,thesis}. Virial coefficients
are calculated in their strong coupling limits. Like the
approximations (equations~(\ref{chixy}) and~(\ref{chiz})) discussed
above, the virial expansion is valid only at sufficiently low
densities that particle chains remain short. When particle chains
exceed a length $L$ equal to the order of the virial expansion
(e.g. $L=3$ in this case) the virial expansion cannot account for
their influence on thermodynamic properties.


\section{Conclusions}

We have developed and tested a cluster-based method for Monte Carlo
simulation of ferrofluids. While we demonstrate this method for
dipolar hard spheres confined in a plane, it should also be applicable
to three dimensions and in the presence of softer short-range
interactions. The method is of greatest utility, compared to
conventional Monte Carlo, when chain formation is prevalent. In the
strong coupling limit our method permits accurate calculation of
magnetic susceptibility that is impossible to achieve in reasonable
time with conventional Monte Carlo or molecular dynamics methods.

Further enhancements of the method may be possible. Note that we
reduced magnetization correlation times more substantially than energy
correlation times. We believe that energy correlation times are
governed by the difficulty of breaking and reconnecting tightly bound
clusters. At present the only means of breaking or reconnecting
clusters is by single-particle motion - we purposely prevent changes
in cluster connectivity during cluster translation due to the
requirements of detailed balance. Conventional Monte Carlo steps
mainly add or remove particles at cluster ends, resulting in only
modest changes in cluster identity. If clusters could be deliberately
broken into shorter segments by multi-particle moves, or joined
together into larger clusters, the energy correlation time could be
further reduced. Because the energy cost of any break is dominated by
a single near-neighbor particle bond, the acceptance rate for such
moves should be about as large as it is for single particle moves.

One possible scheme to incorporate cluster breaking and reconnection
is to use proximity to identify clusters, then execute many cluster
translations during which proximity may change, but always using the
initial cluster identifications when executing translations. Because
cluster breaking is provided subsequent to cluster attachment,
detailed balance can be satisfied.  However, because cluster
identifications are based on past history rather than the
instantaneous configuration, the simulation becomes
non-Markovian. None-the-less, in the limit when cluster identity is
maintained over long time intervals, the simulation should recover the
equilibrium statistics of a genuinely Markov process~\cite{Markov}.

Another possible improvement upon the present method would be to
explore more general types of cluster moves than translation and
magnetization reversal. Inspecting figure~\ref{conf} it is apparent
that undulations of chains are an important type of fluctuation, which
we have not explicitly accelerated using our technique. Such Monte
Carlo steps may indeed be devised~\cite{Sokal} and included in our
simulation, in principle.

In conclusion, our development of a new simulation method enhanced the
accuracy with which magnetic susceptibility may be calculated in the
limit of strong coupling. We applied this simulation technique to
calculate magnetic anisotropies in ferrofluid thin films. Our results
show that a model of magnetic anisotropy based upon clusters confined
within the film suffices to explain the anisotropy provided clusters
are short and well isolated from each other. Both this approximation,
and the virial expansion, fail at high densities when loops and
branched networks form. In this circumstance, computer simulation is
the only method known to calculate magnetic susceptibility. The
cluster-based Monte Carlo simulation described here is needed to
magnetic fluctuations in a reasonable amount of time.

\acknowledgments

This research was supported by NSF grant DMR-9221596 and REU supplements.
We thank H. Zhang and S. Banerjee for useful discussions.

\begin{thebibliography}{99} 

\bibitem{Rosensweig} R.E. Rosensweig, {\it Ferrohydrodynamics}
(Cambridge, 1985)

\bibitem{deGennes} P.G. deGennes and P.A. Pincus, J. Phys. Kondens. Mat.
{\bf 11}, 189 (1970)

\bibitem{smalllambda} P. Debye, Phys. Z. {\bf 13}, 97 (1912);
L. Onsager, J. Am. Chem. Soc. {\bf 58}, 1486 (1936); M.S. Wertheim, J.
Chem. Phys. {\bf 55}, 4291 (1971)

\bibitem{virial} W.H. Keesom, Comm. Phys. Lab. Leiden, Suppl.
{\bf 24 b} (1912); A.D. Buckingham and C.G. Joslin, Mol. Phys. {\bf
40} (1980) 1513

\bibitem{WZ} M. Widom and H. Zhang, in {\it Complex Fluids} eds.
E.D.  Sirota, et. al. MRS conference proceedings {\bf 248}, 235 (1992)

\bibitem{thesis} H. Zhang, {\it Ferrofluids in various forms and phases}
(PhD Thesis, Carnegie Mellon University, 1994)

\bibitem{living} P.C. Jordan, Mol. Phys. {\bf 25}, 961 (1973);
M.E. Cates, J. Phys. France {\bf 49}, 1593 (1988);
J.M. Tavares, M.M. Telo Da Gama and M. A. Osipov, Phys. Rev. E {\bf 56}
(1997) R6252

\bibitem{chainsim} S.C. McGrother and G. Jackson, Phys. Rev. Lett. {\bf 76}
(1996) 4183; M.J. Stevens and G.S. Grest, Phys. Rev. E {\bf 51}
(1995) 5976; M.E. van Leeuwen and B. Smit, Phys. Rev. Lett. {\bf 71}
(1993) 3991; J.J. Weis and D. levesque, Phys. Rev. Lett. {\bf 71}
(1993) 2729; M.J. Blair and G.N. Patey, Phys. Rev. E {\bf 57} (1998) 5682

\bibitem{Ork} G. Orkoulas and A. Panagiotopoulos, J. Chem. Phys.
{\bf 101} (1994) 1452

\bibitem{DHS} For reviews see J.M. Deutch, Ann. Rev. Phys. Chem.
{\bf 24} (1973) 301; M.S. Wertheim, {\it ibid}, {\bf 30} (1979) 471

\bibitem{BC} S.W. de Leeuw, J.W. Perram and E.R. Smith, Ann. Rev. Phys.
Chem. {\bf 37}, 245 (1986); M.P. Allen and D.J. Tildesley, {\it Computer
Simulation of Liquids} (Oxford, 1987).

\bibitem{Weis} J.J. Weis, Mol. Phys. {\bf 93} (1998) 361

\bibitem{Metropolis} N.A. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth,
A.H. Teller and E. Teller, J. Chem. Phys. {\bf 21} (1953) 1087

\bibitem{multipole} L. Greengard and V. Rokhlin, Chem. Scripta {\bf 29A}
(1989) 139; S.W. Yuan and H.N. Bertram, IEEE Trans. Mag. {\bf 28} (1992) 2031

\bibitem{MD} D. Wei and G.N. Patey, Phys. Rev. Lett. {\bf 68} (1992) 2043;
S.W. de Leeuw, C.P. Williams and B. Smit, Mol. Phys. {\bf 65} (1988) 1269;
D. Tomanek, et al. Z. Phys. D {\bf 40} (1997) 539

\bibitem{Markov} D. Bouzida, S. Kumar and R.H. Swendsen, Phys. Rev. A{\bf
45}, 8894 (1992)

\bibitem{Sokal} N. Madras and A.D. Sokal, J. Stat. Phys. {\bf 50}, 109 (1988)

\end{thebibliography} 


\begin{table}
\caption{Temperature dependence of Conventional Monte Carlo correlation times.
All runs with $N=100$ particles at a density of
$\rho^{\star}=0.3$. Correlation times $\tau$ are $1/e$ times.  $L$ is
the mean number of particles per cluster, and $\chi_{xy}$ and $\chi_z$
are the in- and out-of-plane susceptibility, respectively.}
\begin{tabular}{||l|r|r|r|r|r||}
$T^{\star}$ & $\tau_E$ & $\tau_M$ & $L$ & $\chi_{xy}$ & $\chi_z$ \\
\hline
\hline
1.0         & 2        & 5        &    1.5    &  0.5   & 0.20 \\
\hline
0.5         & 8        & 16       &    1.8    &  1.4   & 0.26 \\
\hline
0.3         & 81       & 295      &    3.1    &  3.7   & 0.24 \\
\hline
0.25        & 280      & 3400     &    5.5    &  4.7   & 0.21 \\
\hline
0.2         & 2500     & 110000   &   22.8    &  7.0   & 0.17 \\
\end{tabular}
\label{std_data}
\end{table}

\begin{table}
\caption{Comparison of simulation methods: $C$ denotes conventional; $F$
denotes cluster flips; $T$ denotes cluster translation. $t_{mcs}$ is
time on 50 MHz SPARC 10 in units of CPU seconds per Monte Carlo step.
Simulation temperatures all at $T^{\star}=0.2$. All other parameters
as in table~\ref{std_data}.}
\begin{tabular}{||l|r|r|r||}
method & $\tau_E$ & $\tau_M$ & $t_{mcs}$ \\
\hline
\hline
C      & 2500     & 110000   & 0.22 \\
CF     & 1300     & 100      & 0.65 \\
CFT    & 180      & 53       & 1.26 \\
\end{tabular}
\label{sim_type}
\end{table}

\begin{figure}
\epsfxsize=5in  
\hspace{1in}\epsfbox{conf.ps} 
%\centerline{\epsfig{file=config.eps,width=3in}}
%\vspace{10pt}
\caption{Typical chaining configuration. $N=300$, $T^{\star}= 0.2$ 
and $\rho^{\star} = 0.3$. Because clusters extend across boundaries
of the simulated cell, a 2 $\times$ 2 super-cell is shown.}
\label{conf} 
\end{figure} 


\begin{figure}[tb]
\epsfxsize=5in  
\hspace{1in}\epsfbox{E_plot.ps} 
%\centerline{\epsfig{file=E_plot.eps,width=3.33in}}
%\vspace{10pt}
\caption{Energy autocorrelation functions for
conventional Monte Carlo (solid curve), Monte Carlo with cluster flips
(dashed curve), and Monte Carlo with cluster flips and moves (dotted
curve). $T^{\star}=0.2$ and $\rho^{\star} = 0.3$ as in
figure~\ref{conf}.}
\label{gEofs} 
\end{figure} 


\begin{figure}
\epsfxsize=5in  
\hspace{1in}
\epsfbox{M_plot.ps} 
%\centerline{\epsfig{file=M_plot.eps,width=3.33in}}
%\vspace{10pt}
\caption{Magnetization autocorrelation functions for the same runs as
in figure~\ref{gEofs}. Note that the conventional Monte Carlo method
(solid line) does not decay appreciably over the time scale shown.}
\label{gMofs} 
\end{figure} 


\begin{figure}[tb]
\epsfxsize=5in
\hspace{1in}\epsfbox{anisotropy_plot.ps}
%\centerline{\epsfig{file=anisotropy_plot.eps,width=3.33in}}
%\vspace{10pt}
\caption{Magnetic susceptibility per particle versus density at 
$T^{\star}=0.2$. Data points are simulated, with error bars shown
where appropriate.  Solid and dashed curves are, respectively, second
order and third order virial expansions. Values greater than $\chi_1=1.67$
correspond to $\chi_{xy}$ while values below correspond to $\chi_z$.
Inset: mean chain length $L$ versus density.}
\label{fig:chi}
\end{figure}

\end{document}
