Cluster-based Monte Carlo simulation of ferrofluids

PHYSICAL REVIEW E VOLUME 59, NUMBER 2 FEBRUARY 1999 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, Pennsylvania 15213 ~Received 4 September 1998! 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 two 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. @S1063-651X~99!05702-5# PACS number~s!: 02.70.Lq, 61.20.Ja, 75.50.Mm, 75.70.Ak I. INTRODUCTION Ferrofluids @1# are colloidal suspensions of nanometerscale 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 Fig. 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 Carlo computer simulation, likewise, is difficult when particle chains form @8#. The problem is that the strong attraction of particles within a chain requires 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 multiparticle 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. Mi1063-651X/99/59~2!/2424~5!/$15.00 PRE 59 croscopic 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 m . One defines a reduced temperature T ! 5k B Ta 3 / m 2 and a reduced density r ! 5 r a 3 . The dipole moment of particle i is a three dimensional vector mi , 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 ri 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-two-dimensional 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 DE,0, and accepted with probability exp(2DE/kBT) 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 steps: 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 2424 ©1999 The American Physical Society CLUSTER-BASED MONTE CARLO SIMULATION OF . . . PRE 59 2425 TABLE I. Temperature dependence of conventional Monte Carlo correlation times. All runs with N5100 particles at a density of r ! 50.3. Correlation times t are 1/e times. L is the mean number of particles per cluster, and x xy and x z are the in- and out-of-plane susceptibility, respectively. T! 1.0 0.5 0.3 0.25 0.2 tE tM L x xy xz 2 8 81 280 2500 5 16 295 3400 110 000 1.5 1.8 3.1 5.5 22.8 0.5 1.4 3.7 4.7 7.0 0.20 0.26 0.24 0.21 0.17 angles fall off like AT ! . 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’s!. 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 N5100 particles one conventional MCS takes time t'0.2 sec on a 50 MHz SPARC 10 work station. We find correlation times ~see Table I! in excess of 105 MCS at low temperatures, meaning it requires over 5 h 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 nonphysical 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(2DE/kBT) where DE 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 u ri 2r j u ,a1 e . Typically we choose e 50.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. Fortu- FIG. 1. Typical chaining configuration. N5300, T ! 50.2, and r 50.3. Because clusters extend across boundaries of the simulated cell, a 2 3 2 supercell is shown. ! nately, 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 e 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 pairwise. 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(N21)/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 nonergodic 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 DAVIS, McCAUSLAND, McGAHAGAN, TANAKA, AND WIDOM 2426 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 ! 50.2 and r ! 50.3 as in Fig. 1. 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 ! 5 1 T sim 2t T sim 2t s51 ( @ E ~ s1t ! 2Ē #@ E ~ s ! 2Ē # , ~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 moves ~dotted line!. Figure 3 shows the same for magnetization M (t). Both figures correspond to the density and temperature ( r ! 50.3, T ! 50.2) for which Fig. 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 Figs. 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 of 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 PRE 59 FIG. 3. Magnetization autocorrelation functions for the same runs as in Fig. 2. Note that the conventional Monte Carlo method ~solid line! does not decay appreciably over the time scale shown. such that the typical interchain 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 t 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 x xy and diminishes the out-ofTABLE 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 ! 50.2. All other parameters as in Table I. method C CF CFT tE tM t mcs 2500 1300 180 110000 100 53 0.22 0.65 1.26 CLUSTER-BASED MONTE CARLO SIMULATION OF . . . PRE 59 FIG. 4. Magnetic susceptibility per particle vs density at T ! 50.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 x 1 51.67 correspond to x xy while values below correspond to x z . Inset: mean chain length L vs density. plane susceptibility x z . This behavior is also evident in Fig. 4 which displays the variation of magnetic susceptibility with density r ! at fixed temperature T ! 50.2. We can explain this effect as a result of chain formation. The following discussion presents a semiquantitative analysis. Isolated, noninteracting, particles exhibit Langevin behavior with each particle contributing x 15 m2 3k B T ~2! to the susceptibility in any direction. One factor of m 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 three-dimensional space. For particles in a straight chain of length L, the chain behaves like an effective particle of dipole moment L m and contributes (L m ) 2 /(2k B T) to x xy . The factor of 2 replaces the factor of 3 because the chain is constrained to lie in the xy plane. Collecting N particles into N/L chains of length L the resulting susceptibility per particle is x xy 'L m2 , 2k B T ~3! revealing the enhancement of x xy due to particle chaining. In contrast, the out-of-plane susceptibility x 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 x z' N1 x N 1 ~4! 2427 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 Fig. 4 at low density, below r ! 50.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 Eq. ~3!. Likewise, Eq. ~4! breaks down when so few particles remain unbound that the residual x z contributed by bound particles dominates, and x z becomes independent of density. Above r ! '0.05 the inplane susceptibility x 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 x 1 , given by Eq. ~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 @Eqs. ~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., L53 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 multiparticle 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-nearneighbor 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 2428 DAVIS, McCAUSLAND, McGAHAGAN, TANAKA, AND WIDOM PRE 59 tibility 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. 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. Nonetheless, 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 Fig. 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 suscep- This research was supported by NSF Grant No. DMR9221596 and REU supplements. We thank H. Zhang and S. Banerjee for useful discussions. @1# R.E. Rosensweig, Ferrohydrodynamics ~Cambridge, New York, 1985!. @2# P.G. deGennes and P.A. Pincus, Phys. Kondens. Mater. 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. B 24, Sec. 6 ~1912!; A.D. Buckingham and C.G. Joslin, Mol. Phys. 40, 1513 ~1980!. @5# M. Widom and H. Zhang, in Complex Fluids, edited by E.D. Sirota et al., MRS Symposia Proceedings No. 248 ~Materials Research Society, Pittsburgh, 1992!. @6# H. Zhang, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, 1992. @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, R6252 ~1997!. @8# S.C. McGrother and G. Jackson, Phys. Rev. Lett. 76, 4183 ~1996!; M.J. Stevens and G.S. Grest, Phys. Rev. E 51, 5976 ~1995!; M.E. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 3991 ~1993!; J.J. Weis and D. Levesque, ibid. 71, 2729 ~1993!; M.J. Blair and G.N. Patey, Phys. Rev. E 57, 5682 ~1998!. G. Orkoulas and A. Panagiotopoulos, J. Chem. Phys. 101, 1452 ~1994!. For reviews, see J.M. Deutch, Annu. Rev. Phys. Chem. 24, 301 ~1973!; M.S. Wertheim, ibid. 30, 471 ~1979!. S.W. de Leeuw, J.W. Perram, and E.R. Smith, Annu. Rev. Phys. Chem. 37, 245 ~1986!; M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids ~Oxford, New York, 1987!. J.J. Weis, Mol. Phys. 93, 361 ~1998!. N.A. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 ~1953!. L. Greengard and V. Rokhlin, Chem. Scr. 29A, 139 ~1989!; S.W. Yuan and H.N. Bertram, IEEE Trans. Magn. 28, 2031 ~1992!. D. Wei and G.N. Patey, Phys. Rev. Lett. 68, 2043 ~1992!; S.W. de Leeuw, C.P. Williams, and B. Smit, Mol. Phys. 65, 1269 ~1988!; D. Tomanek et al., Z. Phys. D 40, 539 ~1997!. D. Bouzida, S. Kumar, and R.H. Swendsen, Phys. Rev. A 45, 8894 ~1992!. N. Madras and A.D. Sokal, J. Stat. Phys. 50, 109 ~1988!. ACKNOWLEDGMENTS @9# @10# @11# @12# @13# @14# @15# @16# @17#