Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://chaos.phys.msu.ru/loskutov/PDF/JPA_DLA_2007.pdf Äàòà èçìåíåíèÿ: Mon Jan 31 16:07:54 2011 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 19:56:13 2012 Êîäèðîâêà:

IOP PUBLISHING J. Phys. A: Math. Theor. 40 (2007) 12033­12042

JOURNA

L OF

PHYSICS A: MAT

HEMATICAL AND

THEORETICAL

doi:10.1088/1751-8113/40/40/002

Generalization of the DLA process with different immiscible components by time-scale coarse graining
EBPostnikov1, A B Ryabov2 and A Loskutov2
1 Kursk State University, Theoretical Physics Department, Radishcheva st, 33, 305000, Kursk, Russia 2 MV Lomonosov Moscow State University, Physics Faculty, 119992 Moscow, Russia

E-mail: postnicov@mail.ru, ryabov@polly.phys.msu.ru and loskutov@chaos.phys.msu.ru

Received 28 March 2007, in final form 21 August 2007 Published 18 September 2007 Online at stacks.iop.org/JPhysA/40/12033 Abstract In the framework of the mean-field approximation we propose a new approach to the description of the growth of fractal structures which are formed as a result of the process of diffusion limited aggregation. Our approach is based on the coarse graining of the time scale which takes into account the property of discreteness of such structures. The obtained system of partial differential equations allows us to evaluate numerically the fractal dimension and the cluster density depending on the distance from the cluster center. The results are in a quite good agreement with values found by the direct numerical simulations. The proposed approach is generalized for the case of the cluster description with different immiscible particles. PACS number: 61.43.Hv

1. Introduction The process of diffusion limited aggregation (DLA) is a cluster formation by a sequential adhesion of randomly walking particles. Firstly, this model was suggested in 1981 by Witten and Sander [1], and up to the present it attracts a considerable attention (see, e.g., reviews [2, 3]). This interest is caused by the fact that the DLA process is a suitable model of a number of physical phenomena. For example, the breakdown of dielectric, the growth of viscous fingers or bacterial colonies are similar to a certain DLA process. Moreover, the analysis of DLA clusters as objects with the fractal structure is an interesting mathematical topic. A lot of effort was spent for the numerical investigation of DLA clusters and their scaling properties. In particular, it was found that the density distribution of a DLA cluster has the fractal nature. For example, in 2D case the cluster mass as a function of its radius is a power 1.71 [4]. The process of the diffusion limited growth law M(r) r D with the index D (without long-distance interactions of particles) has been numerically studied, and become even a test platform for the examining new numerical methods [5].
1751-8113/07/4012033+10$30.00 © 2007 IOP Publishing Ltd Printed in the UK 12033


12034

EBPostnikov et al

The great success in the theoretical description of the DLA process was achieved by the fixed-scale transformation, renormalization analysis [6, 7] and iterated conformal transformations [8]. These methods yielded very good analytical estimations of the fractal dimension. However, both methods have the restricted range of applications. The main topic studied via the renormalization analysis is the cluster dimension. But the dynamics of the cluster growth cannot be considered by this method. Iterated conformal transformations are essentially two dimensional. Therefore, a generalization for higher dimensions is not possible. At the same time, in 3D space there exist real DLA clusters: dust clusters and clusters formed by the ion deposition in solutions, etc. Furthermore, some problems of mathematical epidemiology [9] are similar to the DLA process. For example, let us consider an infected individual (a germ of the cluster), and other susceptible individuals, which walk randomly in a vicinity of the germ. For instance, this situation may occur for animals inside a natural habitat as well as for long- or short-range traveling people. After a contact with an infected one, the random walker gets infected and its mobility decreases. Then the distribution of the infected individuals in this model corresponds to a growing DLA cluster. The first mean-field approach [10] for the DLA process was based on relation between the densities of the cluster particles and random walkers. In the ensuing paper [11] a more exact approximation has been proposed. Nevertheless, in spite of the qualitative explanation of the cluster growth, the obtained value of the fractal dimension was far from the real value. That is the reason why in the paper [12] some modifications of the mean-field approach were proposed. For example, it is assumed that the velocity of the cluster growth is proportional to the density gradient of the walking particles. In the other approach [13] an efficient cut-off of the growing cluster front by means of a phenomenological exponent has been introduced. This led to the retardation of the cluster growth. Such retardation allowed to accumulate the particle density in the cluster perimeter and thus, it provided a more accurate value of the fractal dimension. Using an assumption that the cluster and free particles interact with the same probability Bogoyavlenskiy and Chernova [14] showed that the growth velocity is proportional to the second degree of the density and its derivatives. This result can be obtained by the passage from the Boltzmann kinetic integral to a system of differential equations. The particle screening effect and relation of boundary conditions with the cluster surface have been studied in terms of integro-differential equations [15]. It should be also noted that the analysis of the correct boundary conditions has been successfully applied to the description of viscous fingers [16]. Thus, although a certain progress was achieved in the framework of the mean-field approximation and provided a deep insight to the DLA kinetics, it did not lead to the complete solution of the problem. In the present paper we generalize the mean-field approximation in such a way that it makes it possible to obtain an adequate description of the scaling properties of the DLA-cluster growth. Some preliminary results on the fractal pattern formation by DLA process have been described in [17]. Moreover, this approach admits a quite natural generalization for systems containing different immiscible particles.

2. The model In this section we describe the ordinary mean-field approach and present our DLA model.


Generalization of the DLA process with different immiscible components

12035

2.1. Premises of the model In the two-dimensional case the DLA model is formulated as follows. Initially, a fixed `seed' (or germ) particle of the diameter a is located at the center of a circle of the radius R. The second particle appears at the distance R and walks randomly until it hits the seed particle. If the free particle leaves the circle, a new one appears at the distance R from the center. If the particle encounters the seed, then it sticks and becomes a part of the cluster. After aggregation a new particle appears at the binding circle and starts its random walk till the collision with cluster particles, and so on. Thus, during the aggregation process new particles can adhere to any particles of the cluster. To avoid the boundary effects the radius R should sufficiently exceed the cluster size, or R should be increased with the growth of the cluster size. To describe the cluster growth we use the following approximation. Let (r, , t ) be a ^ characteristic function of the cluster, where r and are the polar coordinates and t is time. This function is equal to unity at those points where there are cluster particles, and it is ^ zero otherwise. In the same manner, we define the characteristic function u(r, , t ) of a free particle. In the framework of the mean-field approximation we consider, instead of these functions, the ensemble average. The radial symmetry of the problem allows us to assume that the particle distribution does not depend on the angular coordinate . Thus, ultimately, the cluster density (r, t ) and the density of the free particle u(r, t ) are the functions of only two variables. Additionally, let us define the initial distribution of the free particle u0 (r ) and the zero flux boundary condition ( u/ r |r =0,R = 0). So, the free particle walks randomly until it is aggregated on the cluster surface. From here our approach is a slightly different from [10, 11]. Namely, at each step we look for the probability density of finding a single free particle but not the density of an ensemble of free particles with a constant concentration u(R ) at the binding circle. This approach is more convenient for the time coarsening, which will be introduced later. Moreover, it leads to the same equations in the framework of the ordinary mean-field approximation, with the exception of the outer boundary condition for u(r ). In a totally occupied cell of the same area as a particle the cluster density is = 1/a 2 . However, it is more convenient to take another normalization ( a 2 , u a 2 u), so that or u would be equal to 1 in this case, and characterize the probability of finding a particle within a unit area. Thus, the number of the cluster particles is 1 n = 2 ds, (1) aS where ds is the unit area and S is the area occupied by the cluster. It is more convenient to find the diffusion equation of a free particle in the Cartesian coordinates (x , y ). Let us divide the plane into squared cells of the size a 2 . Then the probability of finding a cluster particle or a free particle in the cell (x , y ) is (x, y, t ) or u(x, y, t ), respectively. Let us consider two independent processes: (1) diffusion of the free particle with its absorption on the cluster; (2) change of the cluster density owing to the absorption of a new particle. Thus, in fact, we neglect the change of the cluster density during the life time of the free particle. Then the difference equation for u(x, y, t ) may be written as follows: u(x - a, y , t ) + u(x + a, y , t ) + u(x , y - a, t ) + u(x , y + a, t ) . u(x, y, t +1) = (1 - (x, y)) 4 (2) If the time step t is small enough then using the Taylor-series expansion of (2), we get u a2 a2 (3) = 2 u - u + 2 u , t 4 4


12036

EBPostnikov et al

where the Laplacian is taken in the Cartesian coordinates. For the radial-symmetric case we have 2 1 2 = 2 + . r r r Up to the terms of the second-order infinitesimal, equation (3) is close to the second equation in the system a2 a2 u a2 (4) = + 2 u = 2 u - + 2 u, t 4 t 4 4 which has been obtained and justified in the paper [11]. One can see a small difference between these equations. On the one hand, in equation (3) the diffusion flux of the free particles is zero in the `solid' cluster area (where n (r ) = 1). On the other hand, if n (r ) = 0, then the last equation of system (4) still allows the aggregation of the free particles on the adjacent sites due to the operator 2 . As was mentioned before, taking into account only the main features of the DLA process, the described mean-field method does not give the correct value of the fractal dimension. The possible reason is that this approach does not consider explicitly the discrete structure of the cluster [18, 19]. 2.2. Time-scale coarse graining In this subsection we propose a mean-field approach for the DLA model based on coarse graining of the time scales. Equation (3) is a diffusion equation with absorption. Thus, the function u(r, t ) decays exponentially with time and the integral U(r ) =
0

u(r, t ) dt

is bounded. This integral, up to a normalization factor Cn , gives the probability density to find a free particle random walking at the distance r until it hits the cluster. Integrating (3) over time from zero to infinity and taking into account that u(r, ) = 0, we arrive at the following equation: a2 2 a2 U(r ) - U(r ) + 2 U(r ) (r), (5) 4 4 where u0 (r ) is the initial distribution of the free particle. The boundary conditions remain the same (i.e. U / r |r =0,R = 0), that means zero particle flux through the boundaries. Now, we turn to the equation for the cluster density (the last one of the system (4)). The change of the cluster density is proportional to the probability of finding the free particle multiplied by the aggregation probability -u0 (r ) = a2 = n+1 (r ) - n (r ) = Cn U(r ) n (r ) + 2 n (r ) , (6) n 4 where Cn is the normalizing factor, and n is the total number of particles in the cluster. In view of (1) and radial symmetry of the cluster, the increment can be written as 2 R (n+1 - n )r dr. a2 0 Form this, the substitution of (6) gives us the value 1 = (n +1) - n = Cn = a 2
R 0 2 a2 4

U(r ) n (r ) +

2 n (r ) r dr

.

(7)


Generalization of the DLA process with different immiscible components

12037

(a)

(b)

Figure 1. (a) The cluster density and the density of equation (8). (b) The cluster density as a function of the (dotted line) and in the model with the coarse graining of the fractal dimension is equal to unit. In the second case

absorbed particles during one step m in radius in the ordinary mean-field model time scale (dashed line). In the first case D = 1.78.

Taking into account the normalizing factor (7), the system of equations (5) and (6) replaces system (4). It is obvious, it has general disadvantages of the mean-field approach. Namely, each new particle slightly increases the cluster radius, that is equivalent to `smearing' of the particle over the cluster perimeter. In the real situation, however, absorbtion of a particle does not change the aggregation probability at the other sites. To construct more realistic equations, let us introduce the following time-scale coarse graining. Let us divide the growing cluster into concentric layers of the width of a particle diameter. We suggest that, the outer layer of the cluster should be completely filled up before the cluster density changes. Afterwards a new outer layer appears, which in turn, should be also fully occupied by the particles, and so on. The number of particles in the layer m of the radius Rm is approximately equal to 2Rm /a . Hence, the normalizing factor Cn should be increased by this number. Unfortunately, we do not know the exact value of Rm . However, let us compare the cluster density (r) and the function U (r )( + 2 ), which shows the aggregation probability (figure 1(a)). The latter function has a clear peak in the neighborhood of the cluster boundary (the border of the outer layer), and equals to zero in the other areas. Thus, in the first approximation, the value Rm can be simply replaced by the variable r. Indeed, due to the peak the magnitude of the right-hand side of (6) vanishes at all r except for r Rm . Then, finally, equation (6) takes the following form: m+1 (r ) - m (r ) = 2r a2 Cn U(r ) m (r ) + 2 m (r ) . a 4 (8)

Now the system consists of three equations. The first one is the ordinary differential equation (5), which should resolve for the cluster of m layers. The second equation describes the normalizing coefficient Cn , which can be found from (7). The recursive formula (8)gives the density of the cluster containing m + 1 layers. Then the procedure recurs. 2.3. Numerical analysis in the plane The numerical solution of equation (5) was based on the finite-difference method using a uniform grid containing 104 points. The initial densities of the cluster and free particles were


12038

EBPostnikov et al

Figure 2. DLA cluster formed by two immiscible particles (light and dark). The inset figure shows initial conditions.

chosen in the form 0 (r ) = exp(-2r/a) and u0 (r ) = exp(-(r - R)2 /a 2 ), respectively. The diffusion terms were discretized by a symmetrical second-order formula. Then, at each step, the obtained tridiagonal system of the algebraic equations was resolved by the backward-sweep method. Integration of expression (7) was made according to the trapezoidal rule. Afterwards, the increment of the cluster density (8) was calculated for the cluster with m + 1 layer. To satisfy the condition of impenetrability of the central part, the increment was taken to be zero at r = 0. Figure 1(b) shows the cluster density obtained from system (5), (8) in comparison with the standard mean-field approximation (4). The latter solution has the slope of about -1, that corresponds to the underestimated value of the fractal dimension (D = 1), whereas our solution (with the slope equal to -0.22 ± 0.02) gives the fractal dimension D = 1.78 ± 0.02. This is very close to the real value. Furthermore, in contrast to the standard approach, which gives a very sharp cut-off of the cluster and the mean-perimeter density drastically falls in zero, our approach leads to the smooth decrease of the cluster density. 3. Aggregation in the system with immiscible components The described mean-field approach with the sequential time coarse graining allows us to generalize the DLA process for system with different immiscible particles. 3.1. Kinetics of the model Consider a 2D area bounded by an impenetrable circle of the radius R.Let N0 be a number of particles of K different kinds initially occupy the center (the inset in figure 2). The arbitrary chosen new particle of a certain kind appears at the border of the circle and walks randomly. If this particle collides with a germ particle of another kind, the random walker will be reflected, and it continues the motion. After the collision with a germ particle of the same kind, they stick together with the probability 1 and a new particle of an arbitrary kind appears at the distance R from the center. We suggest one and the same probabilities pi = 1/K for all kinds of new particles. The described rules lead to the following structure (figure 2). It should be noted that the similar approach has been advanced in the paper [20]. However, the authors considered the diffuse process as a generalization of the problem concerning


Generalization of the DLA process with different immiscible components

12039

percolation in the plane. In their paper, the initial particles were distributed along a line, and the following question was considered. If particles of two kinds are produced with the probabilities p and 1 - p, respectively, what is the value of p providing the infinite growth of `trees' of both kinds. That is the reason why, the growth process was described via a real-space renormalization group analysis. The cluster parts are identified with two possible percolation channels. However, the dynamical properties of the growth have not been considered in detail. Moreover, the fractal properties of the coexisting two families of `trees' remain to be explored. All these questions require a new approach described below. 3.2. Mean-field approach with the sequential coarse graining of the time scale Let us denote the density of free walking particles and the cluster particles as ui (r ) and i (r ), respectively, where i takes the values from 1 to K. A simple modification of equation (5) is necessary to take into account different kinds of particles. Since any particle is a barrier for the other one, the diffusive coefficient should have the same form as in (5). The only difference is a usage of the sum of all densities instead of the single one. But the reaction term must consists of only the functions corresponding to the one kind of particles (i and Ui ). Thus, the system has the following form: K 2 a -ui (r, 0) = j (r ) 2 U i (r ) - U i (r ) i (r ). (9) 1- 4 j =1 Analogously, the cluster growth is described by the system of the independent equations for individual densities:
i i m+1 - m = C m U i i m i m +

a2 2 i m . 4

(10)

i Substitution of m+1 instead of m+1 into equation (7) gives the specific renormalization i constant Cm for each kind of particles, and the factor 2r provides the equal probability for aK the particles of all kinds to occupy the outer layer.

3.3. Numerical simulation: two and three immiscible components This section presents the results of numerical simulation of DLA cluster on a square lattice consisting of two and three immiscible kinds of particles. The initial seed was chosen in the form of an octagon filled with particles along its boundary. The inset in figure 2 shows an example of the initial distribution for the two-kind model. It contains eight separated areas of each type. Our simulation with different initial distributions showed that the final cluster usually has no more than five branches of each kind (figure 2) independently on the initial distribution. Every branch can be divided into two subbranches growing along the branch of the other kind. The fractal dimension of the cluster was obtained by the analysis of the cluster mass M as a function of the radius r (figure 3a). In this figure, the dependence M(r) is shown for the entire cluster and its components consisting of the same kind of particles. Dimensions of the black subcluster and the grey subcluster (see figure 2) are equal to each other (Ds = 1.62). At the same time, the fractal dimension of the entire cluster DK = 1.64 (K = 2). In the case of three immiscible components (see figure 4)wehave D3comp = 1.64 for the complete cluster and Deach = 1.61 for each component. Thus, we may see that the coexistence of different immiscible components leads to the formation of more porous cluster. The cluster


12040

EBPostnikov et al

(a)

(b)

Figure 3. (a) Dependence M(r) for the whole cluster (solid line) and for one component (dashed line). By approximation of linear areas we found fractal dimension to be equal to 1.64 for the whole cluster and 1.62 for one component. (b) Dependences (R ) obtained from equation (9). Fractal dimensions in this case are equal to 1.75 (solid line, K = 2) and 1.78 (dashed line, K = 1).

Figure 4. DLA cluster formed by three immiscible particles.

branches can penetrate to each other. As a result, the dimension of the whole cluster is larger than the dimension of its branches. The result of the numerical solution of equation (9) for K = 2 is shown in figure 3(b) (the black solid line). The gray solid line depicts the linear fit, which allows us to define the fractal dimension. This result (Df = 1.75 for each cluster component) shows that there exists quantitative agreement with the direct numerical simulation. The black dashed line corresponds to the situation with only one kind of particles, and gray dashed line corresponds to its linear fit. One can easily see that the slope of the last line is less, and, respectively, the fractal dimension is more. Analyzing this property, let us compare (5) and (9). It is evident that the diffusion coefficient is sufficiently less in the second case inside the cluster region. For example, in the symmetrical case K = 2 we get 1 - 2 versus 1 - for K = 1. Physically, this corresponds to the constrained motion: an interaction of different particles means the elastic collision, i.e. the particles stop at the given step, and they will be scattered at the next


Generalization of the DLA process with different immiscible components

12041

step. For the identical kinds of the particles this term has the sense of an inelastic collision, i.e. absorption. This absorption will be taken into account by the second term on the right-hand side of (9). 4. Conclusions Thus, in this paper we have presented a continuous model which describes the DLA-fractal growth by a set of difference­differential equations. Similar to the previous investigations (see, e.g., [11, 14, 18]), our approach is based on the radial symmetry of aggregates. However, introducing the coarse graining of the time scale, our approach allows us to take into account the discrete cluster structure. We suppose that the cluster density is discretely changed after the event that the accessible perimeter has been filled completely. As a result, the set of equations describing the aggregate growth has the difference­differential form with the discrete time step. Owing to the fact that in the mean-field approximation the cluster density has an abrupt drop at the boundary (see figure 1), the relation ( + 2 )U has the form which is close to the Dirac -function. This allows us to select the area of the active growth and determine the normalizing factor, which specifies how many particles can be put in the outer layer until it is filled. For two-dimensional case we have a quite good agreement with the known results: D = 1.78 ± 0.02 versus 1.72 in the case of numerical simulations. Moreover, this approach may be easily generalized for the case of larger dimensions. In particular, for a cluster in 3D space we have found that its dimension is about 2.85. Numerical analysis yields D = 2.5. Due to more dense packing of particles this difference is increasing with the growth of the dimension of the embedding space. Furthermore, the developed approach admits a natural generalization for the case of different immiscible components. In this situation we see a good estimation for the cluster dimension (see section 3). The main significance of the proposed approach is not restricted by the description of the DLA process. This method seems to be a quite promising as a tool for the analysis of the contact process, where the result of reaction is determined by the state of agents in neighboring points. Now such processes attract a considerable interest in mathematical epidemiology. However, here integro-differential equations are usually applied (see [9]). At the same time, partial differential equations provide more clear results [21]. The proposed approach has some similar features with the noise-reduced DLA, which was advanced in the paper [22] published in 1986. It has been under active consideration during next years, see the book [23] containing a comprehensive review of the obtained results and approaches. Its main postulate is the following: a random walker is assumed to be stuck to some perimeter point only if the probability of reaching this point exceeds the prescribed level after several attempts. It is found that the relatively small noise-reduction level mainly affects on the branching structure but not on the fractal dimension, if the dimension is defined by scaling of the gyration radius. Furthermore, in the limiting case of zero noise the cluster growth can be reduced to a level-by-level generation [24]. In this approach a number of free particles is added to the growing layer of cluster in the same step. Then the adsorption probabilities are recalculated and a new layer can be filled up. This process is fully determined by the dynamical rules, and the corresponding aggregates have a quite regular form. However, their gyration radii scale like a power function, and the corresponding fractal dimension on the square lattice is D = 1.5 [25]. However, this approach still considers a microscopic level: (i) each site may be filled or not, (ii) only the perimeter sites are considered. At the same time, our proposed


12042

EBPostnikov et al

approach is based on the macroscopic continuous density functions. Thereby, we exclude the explicit stochasticity by averaging over an ensemble of clusters. However, the usage of the differential equations allows us to describe the process of aggregation on the active (outer) layer as well as on layers in front and behind it, corresponding to the aggregation on short and long DLA branches. In addition, the introduction of the level-by-level procedure of the cluster formation gives the fractal dimension which is close to the numerical value obtained for the classical DLA cluster. Finally, we should mention stochastic cellular automata [26]. The correspondence of the description by PDE and cellular automata one can find in [27]. As is known, generation of patterns similar to the DLA clusters has been considered in a pioneer paper [28]. Recently, cellular automata attract a great attention as models of epidemiological processes in inhomogeneous communities (see, e.g., [29]). Our proposed model with the sequential use of a `nonlocal Laplacian' and the secondary averaging seems to be a quite appropriate for the mean-field description of such patterns in radially symmetric case taking into account the cell size (this last condition is a basic one for cellular automata). Acknowledgments The authors would like to thank Professor Nikolai Brilliantov for useful discussions and valuable remarks. References
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] Witten T A Jr and Sander L M 1981 Phys. Rev. Lett. 47 1400 Halsey T 2000 Phys. Today 53 36 Roldughin V I 2003 Russ. Chem. Rev. 72 823 Meakin P 1995 Physica D 86 104 Tillberg D and Machta J 2004 Phys. Rev. E 69 051403 Erzan A, Pietronero L and Vespignani A 1995 Rev. Mod. Phys. 67 545 Hastings M 1997 Phys. Rev. E 55 135 Davidovitch B, Levermann A and Procaccia I 2000 Phys. Rev. E. 62 R5919 Medlock J and Kot M 2003 Math. Biosci. 184 201 Witten T A Jr and Sander L M 1983 Phys. Rev. B 27 5686 Ball R, Nauenberg M and Witten T A Jr 1984 Phys. Rev. A 29 2017 Ohno K, Kikuchi K and Yasunara H 1992 Phys. Rev. A 46 3400 Levine H and Tu Yu 1992 Phys. Rev. A 45 1053 Bogoyavlenskiy V A and Chernova N A 2000 Phys. Rev. E 61 Muthkumar M 1983 Phys. Rev. Lett. 50 839 Feigenbaum M J, Procaccia I and Davidovich B 2001 J. Stat. Phys. 103 973 Ryabov A B, Postnikov E B and Loskutov A 2005 JETP 101 253 Sakaguchi H 1999 J. Phys. Soc. Japan 68 61 Loskutov A, Andrievsky D, Ivanov V, Vasiliev K and Ryabov A 2000 Macromol. Symp. 160 239 ´ Nagatani T and Sagues F 1991 Phys. Rev. A 44 6723 Postnikov E B and Sokolov I M 2007 Math. Biosci. 208 205 Nittmann J and Stanley H 1986 Nature 321 663 Meakin P 1998 Fractals, Scaling and Growth Far from Equilibrium (Cambridge: Cambridge University Press) Batchelor M T and Henry B I 1992 Phys. Rev. A 44 4180 Batchelor M T and Henry B I 1992 Physica A 191 113 Manneville P 1989 Cellular Automata and Modeling of Complex Physical Systems (Berlin: Springer) Loskutov A and Mikhailov A S 1990 Introduction into Synergetics (Moscow, Russian: Nauka) Ulam S 1962 On some mathematical problems connected with patterns of growth of figures Mathematical Problems in the Biological Sciences ed R Bellman (Providence, RI: American Mathematical Society) p 215 [29] Fu S C and Milne G 2004 A flexible automata model for disease simulation Lecture Notes in Computer Science vol 3305 (Berlin: Springer) p 642