Cluster-based Monte Carlo Simulation of Ferrofluids

S.W. Davis, W. McCausland, H.C. McGahagan, C.T. Tanaka and M. Widom

Department of Physics, Carnegie-Mellon University Pittsburgh, Pa. 15213 (October 27, 1998)

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. PACS: 75.50.Mm, 61.20.Ja, 75.70.Ak, 82.70.Dd

Typeset using REVTEX

I. INTRODUCTION

Ferrofluids [1] are colloidal suspensions of nanometer-scale magnetic particles. When the magnetic interactions between particles are stronger than thermal energy, the particles aggregate into chains [2], with the magnetic moment of each particle pointing towards a neighbor along the chain (see figure 1). 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 [3]. Virial expansions [4–6] 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” [5,7], it is important to test such approximate theories with computer simulation.

Traditional Monte Computer simulation, likewise, is difficult when particle chains form [8]. 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 I, 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 [9]. 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 I) due to confinement of clusters within a plane.

II. MODEL AND ALGORITHMS

We model the ferrofluid as a gas of dipolar hard spheres [10]. Each sphere has diameter a and magnetic moment μ. One defines a reduced temperature T^{⋆} = k_{B}T a^{3}/μ^{2} and a reduced density ρ^{⋆} = ρa^{3}. The dipole moment of particle i is a three dimensional vector μ_{i}, and the dipole interaction falls off as 1/r^{3}. To avoid the difficulty of summing the dipole interaction in three-dimensional space [11], we choose to confine our particle positions r_{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 [5]. A related model, in which the dipole moment also is constrained to the plain, has been examined by Weis [12].

We use the Metropolis Monte Carlo method [13], following a Markov Chain through configuration space. Each Monte Carlo step is accepted with probability 1 if the resulting energy change ΔE < 0, and accepted with probability exp (−Δ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^{⋆} << 1. In order for steps and rotations to be accepted with probability near 0.5 (as is the case for the runs in table I) their size must be very small at low temperatures: step sizes fall off proportionally to T^{⋆}, while rotation angles fall off like \sqrt{T^{⋆}}. 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 [14] 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 ≈ 0.2 sec on a 50 MHz SPARC 10 workstation. We find correlation times (see table I) 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 [15], 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 (−ΔE/k_{B}T) where Δ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 |r_{i} − r_{j}| < a + ϵ. Typically we choose ϵ = 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 ϵ 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 2 displays time-correlation functions for energy defined as

g_{E}(t) = \frac{1}{T_{sim} − t} ∑_{s=1}^{T_{sim}−t}(E(s + t) − Ē)(E(s) − Ē)  (1)

where E(t) is the total energy of the configuration after t Monte Carlo steps and Ē 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 3 shows the same for magnetization M(t). Both figures correspond to the density and temperature (ρ^{⋆} = 0.3, T^{⋆} = 0.2) for which figure 1 shows a typical configuration.

Correlation times and run times are tabulated in Table II. 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 2 and 3 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 τ_{M}, accurate calculation of magnetic susceptibility becomes practical, even in the presence of chain formation when conventional methods require absurdly long runs.

III. MAGNETIC ANISOTROPY

With reliable values for the susceptibility in hand, what can we do with the results? Table I reveals growing magnetic anisotropy as temperature drops. Particle clustering enhances the in-plane susceptibility χ_{xy} and diminishes the out-of-plane susceptibility χ_{z}. This behavior is also evident in figure 4 which displays the variation of magnetic susceptibility with density ρ^{⋆} at fixed temperature T^{⋆} = 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

χ_{1} = \frac{μ^{2}}{3k_{B}T}  (2)

to the susceptibility in any direction. One factor of μ 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μ and contributes (Lμ)^{2}/(2k_{B}T) to χ_{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

χ_{xy} ≈ L\frac{μ^{2}}{2k_{B}T}  (3)

revealing the enhancement of χ_{xy} due to particle chaining. In contrast, the out-of-plane susceptibility χ_{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

χ_{z} ≈ \frac{N_{1}}{N} χ_{1}  (4)

where N_{1} is the number of isolated particles not belonging to clusters.

Both approximations (3) and (4) provide excellent approximations to the data in figure 4 at low density, below ρ^{⋆} = 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 (3). Likewise, equation (4) breaks down when so few particles remain unbound that the residual χ_{z} contributed by bound particles dominates, and χ_{z} becomes independent of density. Above ρ^{⋆} ≈ 0.05 the inplane susceptibility χ_{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 [4–6]. The noninteracting Langevin value χ_{1}, given by equation (2), 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 [4–6]. Virial coefficients are calculated in their strong coupling limits. Like the approximations (equations (3) and (4)) 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.

IV. 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 [16].

Another possible improvement upon the present method would be to explore more general types of cluster moves than translation and magnetization reversal. Inspecting figure 1 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 [17] 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.

REFERENCES

[1] R.E. Rosensweig, Ferrohydrodynamics (Cambridge, 1985)

[2] P.G. deGennes and P.A. Pincus, J. Phys. Kondens. Mat. 11, 189 (1970)

[3] P. Debye, Phys. Z. 13, 97 (1912); L. Onsager, J. Am. Chem. Soc. 58, 1486 (1936); M.S. Wertheim, J. Chem. Phys. 55, 4291 (1971)

[4] W.H. Keesom, Comm. Phys. Lab. Leiden, Suppl. 24 b (1912); A.D. Buckingham and C.G. Joslin, Mol. Phys. 40 (1980) 1513

[5] M. Widom and H. Zhang, in Complex Fluids eds. E.D. Sirota, et. al. MRS conference proceedings 248, 235 (1992)

[6] H. Zhang, Ferrofluids in various forms and phases (PhD Thesis, Carnegie Mellon University, 1994)

[7] P.C. Jordan, Mol. Phys. 25, 961 (1973); M.E. Cates, J. Phys. France 49, 1593 (1988); J.M. Tavares, M.M. Telo Da Gama and M. A. Osipov, Phys. Rev. E 56 (1997) R6252

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

[9] G. Orkoulas and A. Panagiotopoulos, J. Chem. Phys. 101 (1994) 1452

[10] For reviews see J.M. Deutch, Ann. Rev. Phys. Chem. 24 (1973) 301; M.S. Wertheim, ibid, 30 (1979) 471

[11] S.W. de Leeuw, J.W. Perram and E.R. Smith, Ann. Rev. Phys. Chem. 37, 245 (1986); M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Oxford, 1987).

[12] J.J. Weis, Mol. Phys. 93 (1998) 361 [13] N.A. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J.

Chem. Phys. 21 (1953) 1087

[14] L. Greengard and V. Rokhlin, Chem. Scripta 29A (1989) 139; S.W. Yuan and H.N. Bertram, IEEE Trans. Mag. 28 (1992) 2031

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

[16] D. Bouzida, S. Kumar and R.H. Swendsen, Phys. Rev. A45, 8894 (1992)

[17] N. Madras and A.D. Sokal, J. Stat. Phys. 50, 109 (1988)

TABLES

TABLE I. Temperature dependence of Conventional Monte Carlo correlation times. All runs with N = 100 particles at a density of ρ^{⋆} = 0.3. Correlation times τ are 1/e times. L is the mean number of particles per cluster, and χ_{xy} and χ_{z} are the in- and out-of-plane susceptibility, respectively.

T^{⋆}  τ_{E}  τ_{M}  L  χ_{xy}  χ_{z}
1.0  2  5  1.5  0.5  0.20
0.5  8  16  1.8  1.4  0.26
0.3  81  295  3.1  3.7  0.24
0.25  280  3400  5.5  4.7  0.21
0.2  2500  110000  22.8  7.0  0.17

TABLE II. 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^{⋆} = 0.2. All other parameters as in table I.

method  τ_{E}  τ_{M}  t_{mcs}
C  2500  110000  0.22
CF  1300  100  0.65
CFT  180  53  1.26

FIGURES

FIG. 1. Typical chaining configuration. N = 300, T^{⋆} = 0.2 and ρ^{⋆} = 0.3. Because clusters extend across boundaries of the simulated cell, a 2 × 2 super-cell is shown.

FIG. 2. 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^{⋆} = 0.2 and ρ^{⋆} = 0.3 as in figure 1.

FIG. 3. Magnetization autocorrelation functions for the same runs as in figure 2. Note that the conventional Monte Carlo method (solid line) does not decay appreciably over the time scale shown.

FIG. 4. Magnetic susceptibility per particle versus density at T^{⋆} = 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 χ_{1} = 1.67 correspond to χ_{xy} while values below correspond to χ_{z}. Inset: mean chain length L versus density.