Particle-cluster aggregation with dipalar interactions

The properties of aggregates generated from an ofF-lattice, two-dimensional, particle-cluster aggregation model with dipolar interparticle interactions have been investigated. The fractal dimension seems to be a monotonically decreasing function of the temperature, between a de6nite value close to 1 at T = 0 and the limit T ~ oo, corresponding to di8usion-limited aggregation of particles with no interaction. Temperature and dipolar interactions are introduced by means of a Metropolis algorithm. By analyzing the orientational correlation function, and what we call the orientation probability density of the direction of the dipoles on the clusters, an ordered state is found at low temperatures. This order diminishes when the temperature increases, due to the disorder induced by the fractal geometry of the aggregates. Our study is extended to the three-dimensional case. A disordered state is found even at T = 0, and a remnant order is shown when analyzing related two-dimensional correlations.


I. INTRGDUCTIQN
The irreversible aggregation process of initially dispersed particles to form fractal clusters (i.e., fractal growth phenomena [1]) has attracted a great deal of interest in the last decade, especially since the development of computer models. There are two basic models of fractal aggregation: particle-cluster aggregation (PCA), the most well-known being the diffusion-limited aggregation (DLA) model of Witten and Sander [2], and clustercluster aggregation (CCA) [4,5]. In particular, numerical simulations performed with CCA models can describe the fractal structure of aerosols and colloids [3], which are in good agreement with experimental results. These computer models constitute coarse-grained simulations of the actual physical processes of aggregation, since they consider, in principle, a very simple short-ranged interparticle interaction: a hard-core, plus an infinite potential well on the surface of the particles, and no interactions for all other distances. For this reason, they are well-suited to simulate aggregation with interactions that strongly decay with distance, but they fail to describe aggregation with long-range interparticle forces.
Some authors have extended CCA and PCA models in order to take interparticle interactions into account. Ansell and Dickinson [6] proposed a CCA model in which particles experience hard-core plus van der Waals interactions by means of a Langevin dynamics algorithm.
Meakin and Muthukumar [7,8] developed a model of reaction-limited CCA [9] with attractive and repulsive power-law interactions. On the other hand, Block et OL [10] presented a model of deterministic PCA, in which particles do not undergo a Brownian motion before being attached to the cluster. Instead, they undergo a deterministic trajectory under a power-law force exerted by the particles already incorporated in the cluster. A particularly interesting case is the aggregation of particles subject to dipolar forces, observed in ferroHuids and the so-called magnetic holes [ll]. Part of its interest lies in the simplicity of the experimental setup for its study, which can be easily performed with magnetized microspheres [12]. This kind of interaction has been simulated in a hierarchical CCA model [13] and in a box with periodic boundary conditions on its sides [14]. In this case, interactions are introduced by means of a Metropolis algorithm [15] and the results obtained are in close agreement with experience [14]. A hierarchical CCA model with ballistic trajectories [16] has also been proposed, yielding results that confirm those in Refs. [13,14]. However, because CCA simulations require a great amount of CPU time to generate large clusters, only modest-sized aggregates can be grown by means of these algorithms (up to 128 particles [13,14,16]). Those results do not allow categorical conclusions to be made about the fractal character of the dipolar clusters obtained.
Our purpose in this paper has been to extend the PCA model in order to include anisotropic dipolar interparticle interactions, by means of a Monte Carlo algorithm based on those described in Refs. [13,14]. Our algorithm is considerably faster than any other CCA dipolar model. At zero temperature, we have been able to generate clusters of up to 5000 particles in a reasonable amount of CPU time. For larger values of temperature, however, the CPU time needed to grow clusters containing more than 1000 particles is exceedingly large. In any case, we have achieved an increase in size of one order of magnitude, allowing an increased accuracy in the results and a reduction of statistical errors over previous studies. We have distributed the paper in the following way. In Sec. II, we describe our computer algorithm for the PCA of dipolar particles in two dimensions. Section III discusses some properties of these dipolar clusters. In particular, we compute the fractal dimension D as a function of temperature and study the order induced on the orientation of dipoles by its reciprocal interactions, which is reBected in two diferent kinds of orientational correlation functions. In Sec. IV

II. COMPUTER. MODEL
In our two-dimensional (2D) off-lattice model, we consider the aggregation process of magnetic particles of diameter d and magnetic moment p, = pu, with p being the magnetic moment strength and m a unit vector oriented along its direction. The dipolar interaction between two particles i and. j, located at the positions t', . and v'~, respectively, is characterized by the energy After a Monte Carlo time step, we compute a new position v' for the random walker. The particle arrives at this position by means of a Brownian motion, which occurs in a direction chosen at random and whose distance is one diameter d. The movement to the new position is performed rigidly, that is, the vector m does not change its orientation during the motion. The energy experienced in the new computed position v' is u» and the total change in the energy due to the movement is AuT --uT -u~. If AuT We start the simulation with a particle, referred to as the "seed, " placed at the origin of coordinates, bearing a randomly oriented tridimensional vector tci, and parallel to the plane of growth. The following particles are released from a randomly chosen position on a circle of radius B;"centered on the seed. The released particles have assigned a vector u; oriented at random. Each particle undergoes a random walk until it either contacts the cluster or moves away from the origin a distance greater than B "q. In this case, the particle is removed and a new one is launched from the circle surrounding the seed. Since magnetic dipole-dipole interactions are long range, the launching radius must be large enough in order to avoid bias in the shape of the cluster due to the initial conditions, i.e. , the position on the launching circle where the particles were released. We have chosen the values B;"=ABe and B "q --AB, where B is the maximum radius of the cluster. The simulations have been carried out with the values A = 2 and e = 5, which are quite reasonable, due to the power-law decay of dipolar interactions. For A = 2, the energy decays on the average one order of magnitude from A (the neighborhood of the cluster) to B;"(the launching distance). A value of A = 4 has also been tested and it provides, within the error bars, the same result for the fractal dimension but the computation time increases dramatically.
The random walks experienced by the incoming particles are affected by the interactions exerted by the particles already attached to the cluster. We have taken this fact into account by using a Metropolis algorithm partially adapted from Refs. [13,14].
Suppose that the cluster is composed by n -1 particles, placed at the points v, , i = 1, . . . , n -1. In some instant of time, the incoming nth random walker occupies the position v and its interaction with the already attached particles is given by the total energy K= d3kgg T (2.5) (K T can be interpreted as some kind of dimensionless temperature, related to the intensity of the interaction. ) In this latter case, the movement is performed with probability p. In the zero temperature limit, K = 0, a movement is accepted if and only if it reduces the dipolar energy. In the high temperature limit, K~oo, all movements are accepted: the particle undergoes a pure Brownian motion and we recover the classic DI A model. After every accepted movement, the moment of the random walker is relaxed, that is to say, it is oriented in the direction of the total field on its position. Indeed this fact assumes that the relaxation time for the orientation of particles is very short in comparison with the movement of the center of mass of the particles. The relaxation phenomenon will turn out to be crucial when analyzing the order of dipoles. The particles stick to the cluster when they overlap one or more particles already incorporated to the aggregate. The new particle is then attached in the point of its last movement where it first contacted one of those particles. After sticking to the cluster, the new attached particle undergoes a last relaxation and its direction is not changed any more.
For K = 0, this algorithm is quite fast. This is due to the fact that dipolar interactions attract the incoming particles towards the cluster, avoiding time-consuming spurious movements far from it. In this case, it is possible to grow clusters up to 5000 particles in a reasonable amount of CPU time. For K & 10, the algorithm becomes very slow, because of the predominance of thermal disorder over magnetic attraction in the Brownian motion. A cluster of 1000 particles at a temperature K~= 10 (the greatest simulated) requires about 8 -10 h of CPU time on an Apollo 9000/730 workstation; a cluster of 1500 particle may take more than a week of CPU time. We have generated and analyzed about 150 clusters containing 1000 particles, for different values of K ranging from 0 to 10; for each value a number between 10 and 40 of aggregates have been simulated. Figure 1 shows some typical clusters of 1000 particles correspond- ing to four difFerent values of K . In order to detect Rnite-size efFects, we have also generated 10 clusters at K = 0 containing 5000 particles, the largest cluster size available with our computer resources.
errors in this paper are at a 95%%uo confid. ence limit). We thus recover, within the error bars, the well-known value D = 1.715 + O.OO4 [17].
The fractal dimension D has been computed for the aggregates grown at each analyzed value of K . Some remarkable results are in order.
(i) K = 0. For 10 clusters of 5000 particles, we get a fractal dimension D = 1.13+0.01, whereas for 40 clusters of 1000 particles, D = 1.20 + 0.02. We can explain this discrepancy by analyzing D as a function of N for the clusters containing 5000 particles, Fig. 2. We notice that D establishes itself on a value around 1.13 for K & 3000, whereas for N & 3000 its value fluctuates and is somewhat larger. At K = 0, a cluster of 1000 particles has not yet reached a stable state and the fractal dimension we compute from it is slightly overestimated. This is also true for the smaller values of K simulated. In any case, the value of D is clearly lower than the one corresponding to true DLA. In Fig. 1(a) we can see a typical cluster with dipolar interactions grown at K = 0. We can compare the form of this aggregate with the one corresponding to true DI A, grown with the same algorithm, which is shown in Fig. 1(d). The dipolar cluster has a much lesser branched and a more open structure, that is to say, it has a smaller bifurcation ratio [18].
(ii) K = 10 . In the cluster-cluster aggregation experiment described in Ref. [14], magnetic particles were employed for which K 1360, at room temperature. In our simulations, the fractal dimension obtained for clusters grown at K i = 10 s is D = 1.35+0.04. Figure 1 shows a typical cluster of this kind. This value of D is clearly difFerent, out of the error bars, from that computed for K = 0.  Figure 3 shows a plot of D as a function of K for the whole range of values analyzed. It seems to exhibit a smooth increasing behavior when increasing K, the values of D ranging between the limits corresponding to K = 0 (aggregation with dipolar forces of infinite strength or zero temperature) and K i = oo (aggregation with no interactions or infinite temperature).
In order to understand this situation we can provide a theoretical reasoning. Let us consider a particle that is going to be attached to the cluster, onto another one belonging to the cluster. It can stick on the tip of a branch with probability Pg, thus contributing to its growth, or it can do so along the branch with probability P, and, therefore, contribute to its splitting.
Both probabilities are related to the dipolar forces and thermal disorder (Brownian motion) through the so-called growth siteprobability distribution [1I. Since dipolar interactions decay as r, we can consider that the main contribution to the energy comes from the interaction between the incoming particle and the particle on the cluster on which the new one is going to be attached. Given the expression (2.2) for the energy, we have two possible positions energetically favored, depicted in Fig. 4, which minimize u;~. Due to the relaxation mechanism, position 4(a) can only happen on the tip of a branch, and therefore it is related to the growth of the branch and to the growing probability P~. Position 4(b), on the other hand, can occur along the whole branch, being related to splitting and to the splitting probability P, . where ug and u, are the energies of the two configurations.
At zero temperature, K -+ oo, we have P, = 0; the most likely event is the growth of a given branch, not its split. Neglecting the effect of thermal disorder, we should then expect at K = 0 an aggregate of D = 1.
In fact, Fig. 2 hints in this direction. D(%) seems to be a decreasing function of N and we cannot reject the possibility of D~1 in the limit N~oo. The actual finite branching found in our simulations at K = 0 is due to the randomness introduced by the Monte Carlo algorithm, which tends to be negligible in the limit of large cluster sizes. For large values of K, we have Pg P" the growth of a given branch and its split are equally likely events. The effect of dipolar interactions is overcome by the Brownian disorder and we recover the original DI A model with no interactions. At intermediate values of K, the ratio between growing and splitting probabilities given by dipolar interactions evolves continuously, as well as the effect of thermal disorder. We should expect a continuous change in the geometry of the cluster and, therefore, a smooth dependence of D on the temperature.

B. Correlation functions
Dipolar interactions, by means of the magnetic relaxation procedure, induce a microscopic order in the orientation of magnetic moments. Neighbor dipoles cannot point in random directions but must be correlated, the orientation of each one being given by the total magnetic field generated by the cluster at the point where the particle was attached. In Fig. 5 we have drawn the orientation of the dipoles associated with the particles belonging to the inner region around the seed, for the clusters shown in Fig. 1. As we can see, dipoles belonging to the same neighborhood point in the direction of the local field on its position, for any value of K One way of obtaining information about the orientation order in the cluster is to study the correlations between the moments of its particles. Following Ref. [19], we define in 2D the ordinary density-density or pair correlation function go (r) = (27r r Ar N) r -~" & I7" -7"'I&r+~~" gg(r) = (2nr&. rN) &" I& + -" where tt(r) is the unit vector oriented along the direction of the magnetic moment of the particle located at position r. This function is a quantitative measure of orientational correlations between pairs of dipoles. We can also define another correlation function of the moments of a cluster g2(r) = (2~rArN) %AT 'MT ) where p(r) is equal to 1 if a particle occupies the position v and 0 otherwise. Lr is the length step in a space discretization. For &actal aggregates, the pair correlation function exhibits a power-law behavior, with gp(r)r the fractal dimension being given by D = d~-a (correlation dimension). In this case, d@ = 2 is the Euclidean dimension of the embedding space.
Following [19], we can define a vector correlation function depending on the relative orientation of pairs of dipoles separated a distance r, I&' -&" I &~+ which depends on the absolute value of the scalar product of pairs of moments. Due to the fractal nature of the clusters, we expect a scaling behavior for these functions [19], gq(r) r ' and g2(r) r~' , and also the same scaling exponent for g2(r) as for the pair correlation vector correlations for dipole orientations show a behavior similar to the correlations for bond vectors described in Ref. [19]. As was found in that paper, gq (r) scales with a different exponent from go(r). However, the difference between the exponents we have found is not as dramatic as it was in the case discussed in [19].
We have also computed the correlation functions from the ensembles of clusters composed of 1000 particles grown at difFerent values of K . From our results we conclude that g2(r) always exhibits the same scaling exponent as the pair correlation function. Figure 7 shows the functions gq(r) computed in a range of values of K between 0 and 10 . These correlations seem to present a power-law behavior with an exponent o. q 1, roughly independent of K . Figure 8 shows gj(r) computed from clusters in the range of K between 10 and 10.
In order to compare the behavior of these functions with a limit case of maximum disorder, we have also plotted in this figure the vector correlations computed from two ensembles of 25 pure DLA clusters (grown with no interactions), each one composed of 5000 particles. The particles of these aggregates have attached dipoles with orientations assigned in two difI'erent ways. For the first ensemble, we assign relaxed dipoles: the seed of the cluster bears a dipole pointing in a randomly chosen direction; the orientation of the particles successively added to the cluster is chosen to be parallel to the total magnetic field in the point where the particle was attached.
For the second ensemble, we assign random dhpoles: the dipole orientation of each particle is chosen completely at random. DLA with relaxed dipoles can be interpreted as the limit case of growth at K~oo with dipolar interactions whose only efI'ect is to relax the vector dipoles in the direction of the total magnetic field. DLA with random dipoles can be seen as the limit case of growth several ensembles of 2D dipolar 1000-particle clusters in a range of values of K between 10 and 10, and from two ensembles of 25 5000-particle 2D DLA clusters. The particles in these latter ensembles have attached relaxed and randomly oriented dipoles, respectively (see text). The vector correlations are short-ranged, with a limit in a maximal disorder situation, DLA with random dipoles. We have used a semilogrithmic plot due to the nonpositive character of gq (r).
with no dipolar interactions, and thus with no ordering efFects at all. In all the cases shown in Fig. 8 Another direct way to measure correlations between pairs of dipoles is simply to study its relative orientation. Consider a pair of dipoles in a cluster whose directions form a certain relative angle 0, which we consider normalized to the interval [0, ir]. We introduce a new quantity, the orientation probability density (OPD) P(K, O), which we define as follows: P(K, 0) dO is the probability that the angle formed by the directions of a pair of dipoles, chosen at random from a cluster generated at some given value This probability density is a measure of the order of the dipoles on the cluster. In a completely disordered distribution, all relative orientations are equally probable and, therefore, we have P(0) = const = -. In a distribution in which all the dipoles point in the same direction, every pair forms a relative angle of zero radians, then P(0) = b(0). A given distribution of dipoles will provide some density P(K, 0) between these two limit values.
In order to determine P(K, 0), we fix a small angular increment 40. On the set of all pairs of dipoles in a cluster (generated at some given value of K and composed of K particles), we compute the function n(K, O), which is de6ned as the number of pairs with a relative angle between 0 and 0+ A9. We then have a'(K) = -) P(K, 0, ). (3.10) In order to reduce statistical uncertainties, these values are averaged over the ensemble of clusters grown at the same temperature. which implies that a(K) = n . A least-squares fitting of P(K, O) on cos0 provides the values for the coefficients a(K) and b(K) given in Table I. As expected, the coefficients a(K) are almost independent of K i and its average over the range of values of K is (a) 3.153 + 0.002, quite similar to the predicted value vr.
For values of K i ) 10  The computed values are given in Table II, and the corresponding average is now (a') = 3.142 + 0.001.
As we can see from our results, P(K, O) can be fitted to the functional form (3.8) for the whole range of K values. The coeKcient a(K) is the normalization constant (a(K) = vr ). The coefficient, which measures the strength of the correlations, is a decreasing function of K; that is, the correlations of relative orientations tend to diminish when increasing the temperature (thermal disorder) towards the limit value b = 0. The limit is already reached for K & 10 The OPD P(K, O) presents quite similar behavior as the vector correlation function gi(r). In the low temperature limit, the orientation of the dipoles is strongly TABLE II. Average value of the OPD of 2D dipolar clusters grown at given values of K 10 1 10 a' (K) 0.3182 + 0.0004 0.3182 + 0.0004 0.3182 + 0.0003 correlated, and this yields a sinusoidal behavior for the orientation density. As the temperature rises, the OPD decreases continuously.
In the high temperature limit, the dipoles show a great deal of disorder (imposed by the disorder of its intrinsic fractal geometry) and the OPD is almost constant (all relative orientations are equally probable). This latter point has been checked by com- Because of limitations in our computer resources, we have only generated clusters of a maximum radius of 75d. This fact imposes an upper limit on the maximum number of particles allowed in a cluster.
In order to check our algorithm, we have grown 10 3D DLA clusters (without dipolar interactions) containing 10000 particles. The radius-of-gyration method yields a fractal dimension D = 2.46 + 0.04, in close agreement within the error bars, with the known value D = 2.495+ 0.005 [17].
We have generated 30 3D dipolar clusters composed of 1000 particles at K = G. Figure 10  The OPD for 3D dipolar clusters P(0) can also be computed by means of expression (3.7). The relative angle 8 in 3D is now the azimuthal angle in spherical coordinates on a relative frame of reference. Figure ll shows P(0) for the ensemble of dipolar clusters, together with the orientation density computed from two ensembles of 3D DLA clusters of 10000 particles, which have assigned relaxed and random dipoles, respectively. All three functions can be fitted to the expression P(g) =sing, which is the density of an equiprobable distribution of azimuthal angles in 3D spherical coordinates. That is to say, in a 3D space the relative orientation of dipoles at zero temperature is randomly distributed. This fact contrasts with the relative order induced by relaxation in 2D. In order to explore this eKect, we have computed the function P'(9), defined as the orientation probability density for the normalized bidimensional projections of the components of the 9D dipoles over fhe X-Y plane. If e stands for an orthonormal base for the X-Y plane, the normalized projection u of the vector n is given by FIG. 12. Orientation probability density for the normalized projections of 3D dip oles over the X-Y plane, P'(8), computed from 3D dipolar clusters and DLA clusters with relaxed and randomly oriented dipoles. The density for 3D dipolar clusters can be 6tted to the function P'(8) = 0.32+ 0.0087cos8 (solid line). We thus recover the sinusoidal behavior shown in 2D, while for the DLA cluster the projected dipoles still present a maximum disorder.
It has no meaning to compute the vector correlation functions of the normalized projections, because they do not have a well-defined minimum cutoff distance (particle diameter in the original 3D space). Figure 12 shows P'(8) computed from the ensemble of 3D dipolar clusters at K = 0, and from two ensembles of 3D DLA clusters with relaxed and randomly oriented dipoles. For the DLA clusters, the OPD of the projected dipoles is almost constant, corresponding to a maximum disorder. However, for dipolar clusters it shows a clear sinusoidal behavior, as it was the case in 2D growth. A least-squares fitting on cos 0 provides the expression P'(8) = 0.32 + 0.0087cos 8. In 3D the dipolar interaction induces some kind of underlying order in the dipole orientations, which only becomes apparent in analyzing the bidimensional properties of the clusters.

V. CONCLUSIONS
In this paper we have investigated the effects of dipolar interparticle interactions in particle-cluster aggregation in two and three dimensions. The fundamental parameter turns out to be the dimensionless magnitude K Eq. (2.5), which is related to the temperature of system T and to the relative strength of the magnetic interactions p. In 2D our results are summarized as follows.
(i) For K i = 0 (low temperature, strong magnetic forces), clusters are formed with a small value of the fractal dimension D = 1.13 + 0.01. The aggregates have a less branched and more open structure than DLA with no interactions. This fact can be explained as the effect of dipolar forces over the growth-site probability distribution of the clusters [1]. The growth of a given branch is much more likely than its split, which results in an enhancement of the screening of the aggregate's inner regions.
(ii) When increasing K (high temperature, weak magnetic forces), we recover a result already found in previous works [13,14,16]: the &actal dimension raises its value until it reaches the limit value of DLA with no interactions. The plot of D as a function of K, Fig. 3, seems to show quite smooth and continuous behavior. This fact can be explained by a continuous relative increment of the splitting probability P, over the growing probability Pg when increasing the temperature, which would cause a rise in the bifurcation ratio and therefore a rise in the fractal dimension. Nevertheless, we cannot utterly reject the possibility of a crossover, with a sharp fall of D at some definite value of K . In order to check this possibility, larger clusters must be generated for intermediate values of K (iii) The short relaxation time for the orientation of dipoles, implemented in our model through the relaxation mechanism discussed, induces a microscopic order in the orientation of neighbor dipoles. At low temperature (K i ( 10 2), this microscopic order induces a macroscopic ordering in the whole cluster, which is revealed by two kinds of measurements: (1) the vector correlation function (3.5) shows a power-law behavior (longrange correlations) with a scaling exponent ni --1.0 roughly independent of K; (2) the orientation probability density (3.7) exhibits a sinusoidal behavior, whose amplitude decreases with K . At high temperatures (K i ) 10 ) macroscopic order is lost, which is shown in the nonpositive definite vector correlations and an almost constant value for the OPD (limit of dipoles randomly oriented). This loss can be explained as the effect of the geometrical disorder induced by the fractal character of clusters. When increasing the temperature, the aggregates are more complex and intricate, and the geometrical disorder prevails over microscopic order, overriding long-range orientational correlations. The ordering properties of dipolar clusters change continuously with K i due to variations in the geometrical disorder (fractal character) of the clusters. This fact suggests a continuous dependence of D on the temperature.
In the 3D case we have only analyzed the limit case of zero temperature. As in the preceding case, a low value for the fractal dimension was found, D = 1.37 + 0.03. In analyzing the orientational correlations, no order is found even at T = 0; in 3D space dipoles are randomly distributed in spite of the relaxation process. Surprisingly, however, when analyzing 2D related properties, namely the OPD of the normalized projections of dipoles over the X-Y plane, we recover the sinusoidal behavior typical of 2D clusters (Fig. 12). In view of these results, we can conclude that, although the relaxation mechanism does not produce an apparent macroscopic ordering in the orientation of 3D dipolar clusters, it indeed induces some underlying order, which is revealed by analyzing the correlations of the normalized projections over a plane. Further work must be done in this direction, including the generation of clusters of size greater than those grown in our present simulations.

ACKNGWI EDGMENTS
The simulations were carried out on the computers of the Centro Europeo de Paralelismo de Barcelona