Furuduy Discuss. Chem. SOC., 1987, 83, 113-124 Computer Simulations of Diff usion-limited Aggregation Processes Paul Meakin Central Research and Development Department, E. I. du Pont de Nemours and Company, Wilmington, Delaware 19898, U.S.A. The diff usion-limited aggregation (DLA) model of Witten and Sander has been used to model a wide variety of physical processes. Here the way in which our picture for the structures generated by DLA models in two dimensions has evolved during the past few years is described. Results from very large scale square-lattice simulations are presented and it is shown how simulations with noise reduction are helping us to understand the effects of anisotropy on the DLA process. It now appears that in the asymptotic (large-mass) limit clusters generated on regular two-dimensional lattices are self-similar fractals with a non-universal fractal dimensionality which is close to but not equal to 1.5. Results are also presented for DLA on two- and three-dimensional percolation clusters.A wide variety of non-equilibrium growth and aggregation phenomena can be described as stochastic processes in which the growth probabilities are determined by a scalar field ( 4 ) which obeys Laplace’s equation (V2@ = 0). Such processes can be simulated using random walks and can often be thought of in terms of particles undergoing Brownian motion. In many important physical realizations the Laplacean field has absorbing boundary conditions ( 4 = 0 ) on the growing structure and a fixed value ( 4 = 1, for example) on a distant surface which encloses the growing objects.Some physical systems which can be described in these terms are given in table 1. The diffusion-limited aggregation (DLA) model was developed by Witten and Sander’ to represent the growth of colloidal aggregates. In this model (illustrated in fig. 1) particles are added, one at a time, to a growing cluster or aggregate of particles via random-walk trajectories. The particles are imagined to be launched from infinity. Since they must cross a circle enclosing the cluster at a random position if growth is to occur, they are in practice launched from a randomly selected point on a circle whose radius is only slightly larger than that of the cluster (fig. 1). In the case of a lattice based model the particles are then transferred to the nearest lattice site and made to follow a random walk on the lattice.Two typical trajectories ( t , and t 2 ) are shown in fig. 1. Trajectory t , eventually moves the particle to an unoccupied surface site (an unoccupied site with one or more occupied nearest neighbours) and growth occurs. Trajectory t2 moves the particle a long distance from the cluster. This trajectory is terminated and a new trajectory is started at a random position on the launching circle to reduce computer time requirements. In fig. 1 the trajectory is terminated at a radius of 3Rm,, where R,,, is the maximum radius of the cluster. In more recent simulations using improved algorithrn~’~~ a termination radius of lOOR,,,,, would be more typical. The termination radius should be sufficiently large so that once a particle has reached the termination circle its position on return to the launching circle will be unbiased, In practice a radius of 3R,,, seems to be sufficient for two-dimensional clusters containing ca.lo4 particles or sites4 and a radius of 100Rm,, is more than sufficient for clusters containing 4 x lo6 sites5 A typical cluster generated using a two-dimensional off -lattice version of this 113114 Diflusion Jim it ed Aggregation Processes Table 1. Some examples of diff usion-limited aggregation process controlling field ref. dielectric breakdown electrodeposition fluid-fiuid displacement in Hele-Shaw cells or porous media materials random dendritic growth dissolution of porous solidification and electric potential or electric potential concentration pressure pressure temperature and/or concentration 29 39-41 19, 21, 33-35 36 37, 38 Fig.1. A simple square-lattice model for diffusion-limited aggregation. This figure shows an early stage in the simulation. At this stage15 (shaded) sites have been occupied including the original seed or growth site (black). Two typical trajectories starting at random positions on the launching circle are shown. Trajectory t , reaches an unoccupied surface site (growth site) and this site is occupied. Trajectory t2 reaches the termination circle which in this case has a radius of 3R,,, where R,,, is the maximum radius of the cluster. This trajectory will be terminated and a trajectory started at a random position on the launching circle. model is shown in fig.2. Considerable interest has developed in DLA and related non-equilibrium growth models, since they lead to the formation of clusters or aggregates having a fractal geometry.' This interest has been sustained by the relevance of the DLA model to a variety of important processes (table 1). The DLA model does not provide a realistic representation of colloidal aggregati~n.~'~ However, it does provide a means of simulating stochastic growth processes controlled by Laplacean (harmonic)P. Meakin 115 4 1000 DIAMETERS b Fig. 2. A 50 000 particle DLA cluster grown using a two-dimensional off-lattice version of the model illustrated in fig. 1. fields. For example, the model illustrated in fig. 1 simulates such a growth process with absorbing boundary conditions on the growing aggregate and a constant field on the termination circle.The growth probability is determined by the local value of the field in the growth sites (or by the field gradient at the surface). Despite the apparent simplicity of the diffusion-limited aggregation model, progress towards a satisfactory theoretical Understanding has been slow. In the absence of such an understanding the interpretation of simulation results is unreliable. However, the results of computer simulations have stimulated and motivated the development of new theoretical ideas9-'* which seem to be leading to a better theoretical understanding of DLA. Our picture for the structure of DLA clusters has undergone considerable evolution during the past few years. The purpose of this paper is to describe some of the recent .advances in our understanding of the structure of DLA aggregates and to indicate some new research directions.Results from Small-scale Simulations The first results from DLA simulations' were obtained using relatively crude algorithms which are quite well represented by fig. 1. Using only slightly improved methods4 it was possible to generate clusters containing up to ca. lo4 particles or occupied lattices in embedding spaces or lattices with (Euclidean) dimensionalities in the range 2-6. These were interpreted in terms of self-similar fractal structures with a charac- teristic fractal dimensionality which was dependent on the dimensionality of the embed- ding space but was insensitive to model details such as sticking probabilities or lattice116 Difusion-limited Aggregation Processes Table 2.Fractal dimensionalities ( 0 ) for a variety of diffusion-limited aggregation modelsa D model 5 0 '/o 75% 2 0 , S(NN) = 1.0 2 0 , S(NNN) = 1.0 2 0 , S(NL) = 1.0 2 0 , S( NN) = 0.25 2 0 , S(NNN) = 0.1 3 0 , S(NN) = 1.0 3 0 , S(NN) = 0.25 3 0 , S(NL) = 1.0 4 0 , S(NN) = 1.0 5 0 , S(NN) = 1.0 6 0 , S(NN) = 1.0 1.73 f 0.06 1.7 1 * 0.08 1.70 f 0.09 1.70 f 0.06 1.70 f 0.10 2.5 1 f 0.06 2.47 f 0.1 5 2.52 f 0.13 3.34f0.10 4.20 f 0.1 1 ca. 5.3 1.70 f 0.06 1.72 f 0.05 1.7 1 * 0.07 1.72 f 0.06 1.69 * 0.08 2.53 rfr 0.06 2.49f0.12 2.50 f 0.08 3.31 f 0.10 4.20 * 0.16 ca. 5.35 structures. At this stage it was generally believed that DLA clusters were homogeneous, self-similar fractals with a broad range of univer~ality.'~ Some of the results used to support this picture are shown in table 2.For two-dimensional clusters using a variety of DLA models (including both lattice and non-lattice models) the radius of gyration (R,) grows with increasing cluster mass ( M ) according to R,- M~ (1) for clusters of small sizes. The exponent appears to have a universal value of ca. 0.585 corresponding to a fractal dimensionality ( Dp = 1/p) of 1.70- 1.7 1 .14 The Effects of Anisotropy The effects of lattice anisotropy first became apparent for two-dimensional square-lattice DLA when improved algorithms273 were developed which allowed clusters containing ca. lo5 sites to be grown. At this size most square-lattice DLA clusters fit into an almost diamond-shaped envelope and it appeared that this might be the asymptotic shape for DLA.This idea motivated the development of a theory for DLA' based on the idea that the asymptotic shape of the cluster is related in a simple way to the structure of the lattice and that the maximum growth probability can be obtained by solving Laplace's equation with simple boundary conditions corresponding to this asymptotic shape. This theory predicted that the fractal dimensionality of DLA clusters should be sensitive to lattice details. While the theory of Turkevich and Scher' is no longer believed to be correct in all of its details, it has motivated much of the subsequent work on DLA, and the idea that the fractal dimensionality of DLA can be obtained from the strength of the leading singularities in the growth probability measure has survived.It still seems probable that the strength of these singularities is related to the shape of the clusters but this shape is not related to the lattice structure in an obvious way. At sizes of lo5 lattice sites or greater it becomes apparent that the overall shape of square-lattice DLA clusters evolves beyond a diamond shape towards a cross-like shape. Fig. 3 shows four maps of the regions in which growth has occurred during the addition of the last 200 000 sites to clusters containing 4x lo6 sites5 This provides a picture of the active zonesI5 (regions where growth is occurring) for the clusters, and clearly illustrates the formation of a cross-like shape. In any event clearer illustration of theP.Meakin 117 h. c -1 I6000 LATTICE UNITS I 6 0 0 0 LATTICE UNITS 4 I6000 LATTICE UNITS 4 > I 6 0 0 0 LATTICE UNITS Fig. 3. Four maps of the active zones for 4 x lo6 site square lattice DLA clusters. Each element in the map consists of a 50 x 50 block of lattice sites and an element is filled if growth has occurred in any of its 2500 lattice sites during the addition of the last 200000 sites to the clusters. effects of lattice anisotropy on square-lattice DLA is provided by the angular distribution of mass (fig. 4). Even for clusters containing lo4 sites the angular distribution is more concentrated along the lattice axes than would be expected for a uniform diamond shape. If clusters are grown on a square lattice with unaxial anisotropy in the sticking probabilities, the cluster has a needle-like shape.The length of the needle in the 'easy' growth directon is given by - M"II (2) and the width of the needle in the 'hard' growth direction (smaller sticking probability) is given by W - Mv-. (3) Ball et a1.l' found that zq = 2/3 and v,. = 1/3. This means that the axial ratio ( Z / w ) diverges as the cluster mass increases and the mean cluster density approaches a constant value (since v,+ vIl = 1). After these results were obtained, it was generally assumed118 Difusion-limited Aggregation Processes Fig. 4. The angular distribution of mass for two-dimensional square-lattice DLA clusters. (a) shows the angular distribution at 5% increments in the cluster mass for clusters containing up to 1.2 x 10“ sites and ( b ) shows similar results for clusters containing up to 4 x lo6 sites.that lattice models for DLA would behave in a similar fashion [i.e. the lengths and width of the cluster arms could be described by eqn (2) and (3)]. A variety of simulation results were found to be consistent with the idea that vIl = 2/3.’4.’6-’8 Results for Y, were much more ambiguous but indicated that 1/3 < v- s 2/3. A value for the exponent vl smaller than vll would indicate that DLA clusters were not self-similar fractals but must instead be described in terms of self-affine fractal geometry.6 The Effects of Noise Reduction A variety of simulation results and have shown that DLA is sensitive to anisotropy and that the shape of DLA clusters is controlled by a ‘competition’ between noise and anisotropy.In the high-noise/low-anisotropy limit the clusters have an irregular shape and the dependence of the radius of gyration on cluster mass is given by eqn (1) with the radius of gyration exponent ( p ) having a universal value of 0.585P. Meakin 119 I I400 LATTICE UNITS I500 LATTICE UNITS I300 LATTICE UNITS I ( d ) 4 I650 LATTICE UNITS Fig. 5. Square-lattice DLA clusters grown with noise reduction parameter (m) of 1 (ordinary DLA), ( b ) 3, (c) 10 and ( d ) 100. for two-dimensional systems. A number of methods have been developed for enhancing anisotropy16-18’22 or reducing noise2716723-25 in DLA simulations. In the most commonly used noise-reduction method the DLA simulations are carried out in the normal fashion except that growth does not occur until a surface site has been reached rn times by the random walkers.After a surface site has been contacted by a random walker the random walk is stopped and a new walker is started from the launching circle after the ‘score’ associated with the contracted site has been incremented by 1. After a surface site has been contacted rn times by random walkers it is filled and any new surface sites are identified and given a score of zero. All of the old surface sites retain their scores, which continue to accumulate as the simulation proceeds. Fig. 5 shows four 50 000 site square-lattice DLA clusters grown with noise reduction parameters (rn) of 1 (ordinary DLA), 3, 10 and 100. The cluster generated with rn = 3120 Diflusion-lim ited Aggregation Processes 750 DIAMETERS < 'CI 1000 DIAMETERS I000 DIAMETERS < + 750 DIAMETERS Fig.6. Results obtained from two-dimensional off-lattice models for diff usion-limited aggregation in which the particles are added in a small number N of directions which are fixed in space. In this figure results are shown for N = ( a ) 4, ( h ) 5, ( c ) 6 and ( d ) 7 with noise reduction parameters ( M ) of 100. shows the effects of lattice anisotropy at least as clearly as the 4~ lo6 site clusters grown with rn = 1. This suggests that the noise-reduction procedure allows us to approach the asymptotic ( M + W) limit for square-lattice DLA without generating enormously large clusters. However, this conjecture has not been firmly established. If this association between large rn and large M is valid then the noise-reduced DLA simulations lead to the following results for square-lattice DLA.( 1 ) In the limit M --+ 00 the clusters are statistically self-similar with a fractal dimensionality which is close to but not equal to 1.5. (2) The axial ratio between the length and the width of the cluster arms ( 1 / w ) approaches a constant limiting value as M + 00. The results have been obtained previously from smaller-scale simulations using the noise-reduced DLA andP. Meakin 121 from related models.16 (3) An angle of ca. 30" can be associated with the tips of the branches of the square-lattice DLA clusters. It appears that it is this angle which determines the maximum growth probability and the fractal dimen~ionality.~~~' An angle of greater than 0" (a finite axial ratio) is consistent with a fractal dimensionality ( D ) of more than 1.5.Similar results have been obtained for growth on a hexagonal lattice with six-fold symmetry and on a hexagonal lattice with growth in three of the six possible directions. In the latter case the results are consistent with vIl = v, (a finite axial ratio in the limit M -+ a) or with vII > vl and vIl - Y, d 0.15. Using lattice models the effects of n-fold anisotropy can be explored where N = 2, 3 , 4 and 6 . A modified off-lattice model (with noise reduction) has been used to explore N = 5 , 7 and 8 also. In this model the vector joining the centre of the contacted particle in the cluster to the centre of the newly arrived particle is-rotated into the closest of N equally separated directions, To implement noise reduction a record is kept of how many times an incoming particle has been rotated to each of the N directions for every particle in the cluster.Fig. 6 shows some of the results obtained for 50000 particle clusters with noise-reduction factors of 100 for the cases N=4, 5 , 6 and 7. For N = 3 , 4 or 5 the clusters always grow with N well defining arms. For N = 6 there are generally fewer than 6 arms using this model, but lattice model simulations with large noise- reduction parameters (m) usually give six well defined arms. It seems that N = 6 is 'marginal' in the sense that the results depend on model details. For N = 7 and 8 results are quite irregular and for N 3 8 they strongly resemble off-lattice DLA clusters. Related Models A variety of models more or less closely related to DLA have been explored during the past few years.This work has been motivated by a need to describe more accurately the events occurring in real processes and by a desire to obtain results which may be used to stimulate and evaluate theoretical work. Much of the present interest in DLA has been the result of its relevance to processes such as fluid-fluid displacement in porous media. Since porous media themselves may have a fractal structure over a significant range of length scales there is reason to investigate DLA on fractal substrates. M e a k i r ~ ~ ~ . ~ ' has simulated DLA on two- and three-dimensional percolation clusters and on a variety of two-dimensional fractals at a finite concentration of random walkers.For two-dimensional percolation clusters at the percolation threshold ( D = 1.89) a fractal dimensionality of 1.40 f 0.05 was obtained, and for three-dimensional percolation clusters ( D = 2.5) a fractal dimensionality of 1.75 f 0.10 was obtained. Murat and Aharony28 have explored the zero concentration limit of this process using both a dielectric breakdown29 version of DLA in which Laplace's equation is solved numerically in the percolation cluster and a model in which the Laplacean field is simulated using random walkers. They find that D = 1.30*0.05 for both two-dimensional models. In a related model3' the growth of the DLA cluster is restricted to a percolation cluster, but the random walkers are not restricted to sites contained in the percolation cluster (they are, however, not allowed to enter sites which are already occupied by the growing aggregate).This process has been simulated in two ways. In the first version of this model the simulation is carried out using a simple modification of the square-lattice DLA model (fig. 1) in which sticking is allowed only at vacant sites which are both nearest neighbours to occupied sites in the growing cluster and are part of an underlying percolation cluster. In the second version the DLA simulation is carried out in the normal fashion except that after a particle (site) has been added to the cluster a random number X uniformly distributed over the range 0 < X < 1 is generated. If X is smaller than a fixed probability, P, the newly added site is sticky and can add more sites.If X a P the newly added site is considered to be 'dead' or 'poisoned' and can add no122 Digusion-Zim i ted Aggregation Processes -500 LATTICE UNITS -& Fig. 7. A cluster of sites grown using the ‘poisoned DLA’ model in which the probability that a site would remain alive (sticky) after addition to the cluster was 0.59. This model is closely related to DLA on a square lattice with growth restricted to a percolation cluster (P, = 0.5927). more sites to its perimeter. For P < P, (the percolation threshold probability) only small clusters grow. As P - P, it becomes possible to grow quite large clusters (fig. 7) and at P = P, this model becomes equivalent to growth on a percolation cluster with unrestricted random walkers. At P = P, a percolation cluster substrate is ‘grown’ as the simulation proceeds by blocking sites with a probability of 1 - P,.If these blocked sites are not included in this cluster the two models become equivalent. This model leads to the formation of clusters which have fractal dimensionality of ca. 1.75. This is not much different than the value of ca. 1.71 found for small square-lattice DLA clusters. However, in other respects the two models are quite different. ( 1 ) Since the underlying percolation cluster is not anisotropic, the clusters are not sensitive to the anisotropy of the square lattice. (2) For square-lattice DLA the dimensionality Dmin, which indicates how the minimum path measured on the cluster between two points on the cluster scales with their separation r measured in the Euclidean embedding space, has a value of 1:31 (4) 1 - yDmln* For DLA with growth restricted to a percolation cluster Dmin must be at least as large as D m i n for a percolation cluster; a value which is believed to be greater than l.32 (3) In contrast to ordinary square-lattice DLA the seed or growth site is not necessarily near the centre of the cluster.If the random walkers are allowed to enter sites which are already occupied by the growing aggregate (from unoccupied sites which are not part of the percolation cluster) then the percolation cluster becomes uniformly filled ( D = 1.89). In this case the trapping of random walkers at the sites on the percolation clusters which are adjacent to filled sites on the aggregate is not sufficient to prevent random walkers from penetrating into the interior of the aggregate.Summary The evolution of square-lattice DLA clusters from a more or less irregular circular shape to a cross-like shape results in different effective exponents (v, and vII) which describeP. Meakin 123 the growth of the width and length of the cluster arms. However, simulations carried out using modified DLA algorithms in which the effects of noise are reduced indicate that the axial ratio (Z/w) reaches a constant value and consequently that there is no divergence of length scales (I/(( = vl in the M -B 00 limit). Consequently, it now seems that the asymptotic geometry for DLA is self similar. Since DLA is relevant to a wide variety of physical processes a large variety of models more or less closely related to DLA but designed to better represent real systems have been developed.Results for DLA on percolation clusters are presented. In this case the substrate has no anisotropy and the clusters grown on percolation clusters (with random walkers both restacked to the percolation cluster and unrestricted) give clusters which are different from ordinary DLA in several important respects. The work described here would not have been possible without contributions from many colleagues. I would particularly like to thank R. C. Ball for his collaboration in much of the work presented here. References 1 T. A. Witten and L. M. Sander, Phys. Res. Lett., 1981, 47, 1400. 2 R. C. Ball and R. M. Brady, J. Phys. A , 1985, 18, L809. 3 P. Meakin, J.Phys. A , 1985, 18, L661. 4 P. Meakin, J. Phys. A, 1983, 27, 604; 1495. 5 P. Meakin, R. C. Ball, P. Ramanlal and L. M. Sander, Phys. Rev. A, submitted for publication. 6 B. B. Mandelbrot, The Fractal Geometry of Nature (W. H. Freeman, San Francisco, 1982). 7 P. Meakin, Phys. Rev. Lett., 1983, 51, 1119. 8 M. Kolb, R. Botet and R. Jullien, Phys. Rev. Lett., 1983, 51, 1123. 9 L. Turkevich and H. Scher, Phys. Rev. Lett., 1985, 55, 1026. 10 R. C. Ball, R. M. Brady, G. Rossi and B. R. Thompson, Phys. Rev. Lett., 1985, 55, 1406. 11 T. C. Halsey, P. Meakin and I. Procaccia, Phys. Rev. Lett., 1986, 56, 854. 12 R. C Ball, in Statistical Physics, ed. H. E. Stanley (North Holland, Amsterdam, 1987). 13 T. A. Witten, in Workshop on Dynamics of Macromolecules, Institute for Theoretical Physics, University of California at Santa Barbara, J.Polym. Sci., Polym. Symp. 73, ed. S. F. Edwards and P. A. Pincus, p. 7. 14 P. Meakin Phys. Rev. A, 1986, 33, 3371. 15 M. Plischke and Z. Racz, Phys. Rev. Lett., 1984, 53, 415. 16 J. Nittmann and H. E. Stanley, Nature (London), 321, 663. 17 J. D. Chen and D. Wilkinson, Phys. Rev. Lett., 1986, 55, 1892. 18 P. Meakin, Z-Y. Chen and P. Evesque, unpublished work. 19 E. Ben-Jacob, R. Godbey, N. D. Goldenfeld, J. Koplik, H. Levine, T. Mueller and L. M. Sander, Phys. 20 M. Matsushita and H. Honda, J. Phys. SOC. Jpn, 1986, 55, 2483. 21 J. Feder, T. Jossang, K. Jorgen and U. Oxaal, in Fragmentation, Form and Flow in Fractural Media, 22 J-D. Chen and D. Wilkinson, Phys. Rev. Lett., 1985, 55, 1892. 23 C. Tang, Phys. Rev. A, 1985, 31, 1977. 24 J. Kertesz and T. Vicsek, .I. Phys. A, 1986, 19, L257. 25 B. R. Thompson, preprint. 26 P. Meakin, in Kinetics ofAggregation and Gelation, ed. F. Family and D. P. Landau (North Holland, 27 P. Meakin, Phys. Rev. B, 1984, 29, 4327. 28 M. Murat and A. Aharony, Phys. Rev. Lett., 1986, 57, 1875. 29 L. Niemeyer, L. Pietronero and H. J. Wiesmann, Phys. Rev. Lett., 1984, 52, 1033. 30 P. Meakin and Y. Termonia, unpublished. 31 P. Meakin, I. Majid, S. Havlin and H. E. Stanley, J. Phys. A , 1984, 17, L,975. 32 S. Havlin and R. Nossal, J. Phys. A , 1984, 17, L427. 33 L. Patterson, Phys. Rev. Lett., 1984, 52, 1621. 34 J. Nittmann, G. Daccord and €3. E. Stanley, Nature (London), 314, 141. 35 H. Van Damme, F. Obrecht, P. Levitz, L. Gatineau and C. Laroche, Nature (London), 1986,320, 731. 36 G. Daccord and R. Lenormand, Nature (London), in press. Rev. Lett., 1985, 55, 1315. Ann. Isr. Phys. SOC., 1986, 8, 531. Amsterdam, 1984), p. 91.124 Difusion -1im ited Aggregation Processes 37 W. T. Elam, S. A. Wolf, J. Sprague, D. U. Gubser, D. Van Vechten, G. L. Barz Jr and P. Meakin, 38 H. Honjo, S. Ohta and M. Matsushita, J. Phys. Soc. Jpn, 1986, 55, 2487. 39 M. Matsushita, M. Sano, Y . Hayakawa, H. Honjo and Y. Sawada, Phys. Rev. Lett., 1984, 53, 286. 40 M. Matsushita, Y . Hayakawa and Y. Sawada, Phys. Rev. A, 1985, 32, 3814. 41 R. M. Brady and R. C. Ball, Nature (London), 1984, 309, 225. Phys. Rev. Lett., 1985, 54, 701. Received 8th December, 1986