Anisotropic percolation and the d-dimensional surface roughening problem

Physica A 200 (1993) 200-211 North-Holland SDI: 0378-4371(93)E0224-3 Anisotropic percolation and the d-dimensional surface roughening problem Sergey V. Buldyrev”, Shlomo Havlin”‘b and H. Eugene Stanley” “Center For Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA bDepartment of Physics, Bar-llan University, Ramat-Gan 52900, Israel We review recent numerical simulations of several models of interface growth in d- dimensional media with quenched disorder. These models belong to the universality class of anisotropic diode-resistor percolation networks. The values of the roughness exponent (Y= 0.63 ? 0.01 (d = 1 + 1) and (Y= 0.48 2 0.02 (d = 2 + 1) are in good agreement with our recent experiments. The values of (Y in higher dimensions (a = 0.38 5 0.03 in d = 4 and [Y= 0.27 ? 0.05 in d = 5) do not support a recent theoretical conjecture. 1. Introduction The growth of rough interfaces in random media is a topic of current interdisciplinary interest [l-5]. For the most part, two types of models have been applied to interface-roughening phenomena: (a) nonlinear Langevin equations - principally the KPZ equation [6] - resulting in self-affine interfaces, and (b) spreading and invasion percolation models [7], generating self-similar interfaces. The self-affine interface (a) can be characterized by the rms surface width w(e, t) = ([h(x, t) - (h(x, t))12)1’2. (1) Here h(x, t) is the surface height at time t, and the angular brackets denote the average over x belonging to a (d - 1)-dimensional hypercube of size gd-’ in the horizontal cross section perpendicular to the direction of growth. An alternative and equivalent quantity is the height-height correlation function ~(8, t). The scaling exponents obtained from w and c are believed to be identical, so we use them interchangeably. Analysis of the KPZ equation implies the scaling law [l-4] w(t, t>- 46 t) - eaf (S), Elsevier Science Publishers B.V.  S.V. Buldyrev et al. I Anisotropic percolation and surface roughening 201 where (2b) and z=ldp. (24 The roughness exponent is cy = l/2 and the dynamical exponent is p = l/3 for d = 1 + 1 [6]. N umerical studies of the KPZ equation give (Y= 2/(d + 2) in d 3 2 [g-lo]. On the other hand, approach (b) (percolation-type models) produces self- similar (fractal) interfaces [7], for which ~(4, m) - f? with (Y= 1 and p = 1. For many phenomena in d = 1 + 1 - from bacterial growth [ 111 and viscous flows [12,13] to the wetting [13-E] and burning [16] of paper-self-affine surfaces are found with anomalous exponents (Yand /3 significantly larger than the KPZ values but less than 1. Recent experimental data in d = 2 + 1 also show anomalously large values of (Y-for mountain surfaces (Y= 0.58 [17,18], for wetting of porous media (Y= 0.5 [19], and for ion beam erosion of metal surfaces (Y= 0.53 [20]. For the past two years, the following question has been addressed: Do these experiments, in fact, represent a crossover from self-affine to self-similar behavior, or is there a new universality class of growth models that produces self-affine interfaces with an anomalous (Y> 2/(d + 2)? As a step in answering this question, we develop and generalize several models [13,14,21-241 of spreading percolation with anisotropy in the growth direction which belong to a new universality class characterized by a self-affine interface with the roughness exponents given in table I. Together with J. Kertesz [25], we have been exploring the relation of this universality class to Table I Critical exponents and percolation thresholds pr, p:, p:, for models A, B and C. d 1+1 2+1 3+1 4+1 0.63 + 0.01 0.48 k 0.03 0.38 * 0.4 0.27 2 0.05 0.63 k 0.01 0.41 k 0.02 0.28 k 0.3 0.18 k 0.03 I 1.01 ? 0.02 1.16 k 0.03 1.32~ 0.5 1.50~0.10 US”,” 1.46 5 0.02 2.18 t 0.03 2.54 2 0.05 3.00 2 0.20 IS 0.60 5 0.03 1.14 f 0.06 1.6kO.l 1.920.2 P: 0.4698 -c 0.0002 0.7423 r 0.0002 0.8423 2 0.0002 0.889 -c 0.0005 PCB 0.5388 2 0.0002 0.8009 t 0.0002 0.8857 2 0.0004 0.9240? 0.0005 c PC 0.6447~ 0.0001 0.9340 + 0.0003 0.9863 2 0.0004 0.9970 k 0.0005 202 S.V. Buldyrev et al. I Anisotropic percolation and surface roughening the reversed percolation transition point [22-241 in the diode-resistor percola- tion network, at which the backflow current emerges in the direction opposite to the direction of diodes, when the concentration of resistors q approaches some critical value q, = 1 -p, (p = 1 - q is the concentration of diodes). In d = 1 + 1, this transition point is dual to the directed percolation point [22-241 and the roughness exponent (Yof the interface between blocked and unblocked diodes can be expressed in terms of correlation exponents of directed percolation. In d > 2, little is known about the critical properties of this transition. Recently, several theoretical conjectures have appeared concerning the values of the critical exponents of the surface roughening models with depinning transition caused by the quenched noise for continuum systems [26,27]. The numerical data for the continuum models in d = 1 + 1 strongly support the argument that these models belong to the same universality class as discrete models described below. However, according to [26,27] d = 5 should be the upper critical dimension above which a = p = 0. Our numerical results do not support either of these conjectures. 2. Models The models A, B and C we have studied are straightforward generalizations to d dimensions of the d = 1 + 1 models defined in [13,14,21,22]. In each model the porous media are simulated by a cubic lattice with certain fraction p of randomly blocked cells which represent the inhomogeneities of the media (see fig. 1). The horizontal cross section of the lattice is the (d - 1)-dimensional hypercube of volume L*-’ with periodic boundaries. Every lattice cell can be wet or dry. At t = 0, all lattice cells are dry except those with vertical coordinate h G 0. At every time step, we simultaneously examine all dry cells on the wet-dry interface and decide whether each of these cells should become wet on the next time step. The decision concerning each cell is taken according to the deterministic rules specific to each model. In analogy with spreading percolation [2,29], we define the tth shell as the set of cells that become wet on the tth time step. In model A, each cell adjacent to the interface becomes wet if (i) it is unblocked or if (ii) it lies below the unblocked cell, adjacent to the interface. In model B, the cell becomes wet if (i) it is unblocked and its nearest neighbor from below is wet or (ii) if the height of the highest wet cell in one of the nearest neighboring columns is larger then the height of the cell under consideration. In model C, rule (i) is the same as in model B, but rule (ii) is slightly  S.V. Buldyrev et al. I Anisotropic percolation and surface roughening 203 Fig. 1. Explanation of the surface growth models A, B and C in d = 1 + 1. Cells are randomly blocked with probability p (indicated by shaded area) or unblocked with probability 1 -p (indicated by a white area). Wet cells are marked by numbers indicating time step at which they have become wet. The line connecting the centers of the dry blocked cells indicates the spanning path of directed percolation that has stopped the growth. (a) Model A, proposed by Buldyrev et al. [13,14] (p, = 0.4698); (b) model B proposed by Tang and Leschhorn [21] (p, = 0.5390); (c) model C proposed by Dhar et al. [22] (p, = 0.6447). The configuration of blocked cells is the same for all three models, thus the largest cluster of model C corresponds to the largest p,. In model C diodes correspond to blocked cells resistors to the unblocked cells. The order of wetting differs from the original definition of ref. [22] in order to prevent overhangs in the moving parts of the interface. These changes. however, do not affect the shape of the pinned interface, which in both cases is the path of the bond directed percolation of diodes on the dual lattice. different. The blocked cell with the coordinates (x1, x2, . . . , xd_*, h) becomes wet if its nearest neighbor from below is wet and at least one of the cells with the coordinates (x1 + Ax,, x2 + Ax,, . . . , .x~_~ + Ax~_~, h + Ah) is wet, where the increments Axi and Ah obey the following constraints: Ah 2 0, lAxi G 1; if Ah >O then it is sufficient that lAxi\ =S1; if Ah = 0 then Axi = 0 or Axi = (-l)h+*Y Model C is in fact equivalent to the diode-resistor network in which each blocked cell corresponds to the diode, and each unblocked cell to the resistor connecting a pair of the opposite vertices of the cell. Each diode or 204 S.V. Buldyrev et al. I Anisotropic percolation and surface roughening resistor has orientation (AX,, Ax,, . . . , Axd_ 1, -l), where Axi = -( - l)‘~+~ are determined by the coordinates of the cell (see fig. lc). Model A is analogous to spreading percolation (rule (i)) with erosion of overhangs (rule (ii)). In models B and C, the height of each column at each time step can increase only by one cell. In model A, the erosion of overhangs corresponds to much faster local growth. Also, in models B and C, the maximal difference in heights of the two neighboring columns is 2, while in model A, it can be any number. Model A was suggested by Buldyrev et al. [13-151, while model B was originally proposed by Tang and Leschhorn 1211 and can be considered as the discretized version of the Langevin equation with quenched noise and is very close to the models defined by Parisi [25] and Csahok et al. [26]. Comparison of the time development of all three models is shown in fig. 1 for the same d = 1 + 1 configuration of blocked cells. The growth of each model is stopped by the spanning paths of directed percolation of the blocked cells with different definitions of connectivity (shown as a fence connecting the centers of the blocked cells). In model A the path has five choices (North, South, East, North-East and South-East), in model B three choices (East, North-East and South-East), in model C only two choices (North-East and South-East). This results in different percolation thresholds: pp = 0.4698, p: = 0.539, p: = 0.6447. Numbers on the cells show the time at which they become wet. The finite size corrections to scaling for these models have different values, but the asymptotic behavior is characterized by the same set of exponents. 3. Dimension 1 + 1: theory and simulations In this section, we review theory and simulations for d = 1 + 1 [13-151. When the probability of blocked cells p is close to p,, the growth is halted in many places by the paths of a directed percolation cluster. Each path can be characterized by two correlation lengths, 5, and E,,. When the path is spanning, i.e. when the growth is stopped completely, t,, is equal to the system size L, and 5, is proportional to the width of the interface w. It is known from the theory of directed percolation [28,29] that the correlation lengths diverge in the vicinity of p,, 5rIP-PcI-y~ Y 5,, - IP -PclrY” ) where y, = 1.733, vI = 1.097 [29]. Thus, W)  S.V. Buldyrev et al. I Anisotropic percolation and surface roughening 205 where (Y= V~/y, = 0.633 + 0.001 . Thus the spanning path of directed percolation describes the final state of the paper wetting experiment [13,14], where the wetting front is completely pinned by inhomogeneities in the paper. The theoretical value of (Y (eq. (3~)) is in excellent agreement with both our simulations and our experiments [13-151. A natural question is “What is the dynamics of wetting?” To answer this, we study the dynamical behavior of the models below and above p,. Fig. 2 shows a snapshot of the wetting front as it continues to propagate in the (1 + l)- dimensional media when p <p,. Large sections of the interface are already pinned and the growth is occurring only in columns that contain unblocked cells on the wet boundary (shown as dark vertical lines). The average horizontal size of pinned sections is e,, , while the average vertical size of pinned sections is 5,. The moving parts have constant steep slope (in model B it is exactly equal to 2), and their vertical and horizontal dimensions are propor- tional to 5,. As was shown in [25], the height-height correlation function c(f?, t), calculated for such an interface for 4 -5, is proportional to e(3u-1)‘2a = f”.71. The effective roughness exponent of this moving interface calculated directly from our numerical data near criticality in d = 1 + 1 is (Y+ = 0.70 ? 0.05, which is in good agreement with the above geometrical arguments. The growth is now mostly a fast erosion of steep slopes, propagat- Fig. 2. Two successive snapshots of the interface of model A still evolving near its pinning threshold (p = 0.4698). System size is 212. (But the calculation is made for L ~2l’.) Upper snapshot (a) taken at t = 18000, lower (b) at t = 19500. Light color indicates wet area, dark color indicates dry area, dark strips indicates “live” columns, i.e. columns that contain cells which become wet at the current time step. Arrows show the direction of the propagation of steep eroding slopes, which has approximately equal horizontal and vertical dimensions both of the order of 5,. The horizontal dimension of the long blocked region is of the order of [,, , while its vertical dimension is of the order of 5,. 206 S.V. Buldyrev et al. I Anisotropic percolation and surface roughening ing horizontally with constant speed. This observation implies that the dynamical exponent zdyn = cq+, l/3 has a value close to 1. Thus /3 = CQ~,,,in good agreement with our numerical results [16-181. This large value of Q,,, may explain the large values 0.7-0.8 found in dynamical experiments [ll-161. In the case d = 2 + 1 we find that for all studied models (Y= 0.48 + 0.02, /3 = 0.41? 0.02, z = 1.16 ? 0.02, while the results for p, are different (see table I). In d > 2 the duality to directed percolation breaks down and the growth of the surface is pinned by the self-affine hypersurfaces, which can be character- ized by horizontal and vertical correlation lengths t,, and 5,. (See ref. [30] for an isotropic model of percolating hypersurfaces.) Below criticality moving parts form steep circular terraces, surrounding pinned parts. However, the moving parts in d > 1 + 1 do not move along straight lines, but rather perform a kind of correlated random walk, which yields z > 1. At d = d,, z should become 2, as for uncorrelated random walk. So far, we find that even for d = 5, z = 1.5 + 0.1~ 2, suggesting that d, > 5. The effective roughness exponent found for the moving interface in d = 2 + 1 near criticality is about 0.52 + 0.05, which is in good agreement with our experiments on wetting of 3D porous media [19]. 4. Avalanches and fractal dust Above p,, the growth is stopped by the spanning path of a directed percolation cluster in d = 1 + 1, or by a self-affine surface in d = 2 + 1. However, we can modify our models and assume that even when the growth is completely stopped, the blocked cells on the interface may still erode-but at an infinitesimal rate. With this assumption, we can remove blocked cells at random when the interface is completely stopped. Each removal will produce an avalanche of growth which eventually will die out when front reaches another directed spanning path, or directed surface of blocked cells (see fig. 3). The distribution of avalanche sizes P(V) is found to be [15] P(V) - v-T~=‘F(v/vo) ) (44 where V is the number of sites removed in an avalanche, and V, - L$‘~, is the characteristic volume. The probability P(V) is estimated to be the ratio of the number of avalanches of size V to the total number of avalanches. In d = 1 + 1, the maximum linear extent of the avalanches (fig. 3), in the longitudinal and transverse directions, is found to scale with exponents zqal= 1.73 * 0.02 ) vyal = 1.10 -+0.02 ) (4b)  XV. Buldyrev et al. I Anisotropic percolation and surface roughening 207 Fig. 3. Successive series of pinned interfaces of model A, showing the boundaries of avalanches, produced by removing a randomly chosen blocked cell from the previously pinned interface. L = 400, p = 0.5 >p,. Correlation lengths r,, and 5, are the typical sizes of the avalanches. in excellent agreement with the correlation-length exponents of directed percolation. Moreover, we find 7 aval= 1.245 + 0.02. (4c) From avalanche studies in d = 2 + 1 we find the correlation-length and roughness exponents to be v,, = 1.06 + 0.1, vI = 0.47 + 0.1, which gives slightly smaller values of LYthan from analysis of height-height correlation function and width: (Y= V, /v,, = 0.44 + 0.1. This may be due to large error bars that are caused by comparatively small errors in the value of p,. An alternative way of producing avalanches is to start growth from a single unblocked cell at time t = 0 when interface is flat (see fig. 4). Above p,, the clusters of wet cells will be all pinned by the blocked cells. Below p, some of these clusters will grow infinitely, but some will be stopped by the pinning surfaces. In analogy with conventional percolation, the survival probability P,,,,(t) of the clusters for t < t, will decay as a power law, P,,,Jt) - t’-7surv. Here t, is the characteristic time which is related to the correlation lengths: 511-G, 5, -t{. The exponent 7S,rVis related to the T,,,, of eq. (4): (r,,,, - 1) = (T,,,, - l)(d - 1 + a)/~ . (54 For t > to,J’,,,, (t) either goes to zero exponentially for p >p, or approaches a 208 S.V. Buldyrev et al. I Anisotropic percolation and surface roughening Fig. 4. Horizontal projection of the cluster of model B (p = 0.8009 =p,) which was started at the center of the screen 21° time steps ago. The current diameter of the cluster is about 21°. The blue area shows the flat interface that is left dry since the beginning of the process. Darkest shades of gray correspond to the largest heights of the interface. Red dots forming “fractal dust” indicate cells that become wet at the current time step. constant value P(m), the probability of an infinite cluster, for p <p,. Thus, studying P,,,,,(t) provides very accurate method of estimating p,. The critical exponents can be derived from the parameters of the clusters exactly at p,. One of these exponents is the exponent S that characterizes the time  S.V. Buldyrev et al. I Anisotropic percolation and surface roughening 209 dependence of the size of the percolation shell: n(t) - t’. It is possible to connect this exponent with z and (Y: s + 1= (d - 1-t a)/2 . (5b) In d = 1 + 1, since z = 1, 6 = (Y= /?, which agrees with the simple geometrical picture (see fig. 2) that the projection of the shell corresponds to the length of the largest steep moving terrace, which scales as the vertical size of the whole system w(t) - to. We have studied the time dependence of the size of the shell n(t) - t* numerically. The results are in agreement with the above relations. The projection of the shell forms a fructul dust (see figs. 2 and 4). In d = 1 + 1, when z = 1, the fractal dimension of this dust d, is connected to the exponent of distributions of avalanches (i.e. blocked regions) through the relation d, = rsurv - 1 = 0.46 ? 0.02. (6) Note that d,< 6, which means that, in fact, fractal dust is packed in moving blocks (the largest moving block is about of the vertical system size w(t)). These moving blocks behave like quasi-particles which are distributed in a fractal way with fractal dimension equal to d,. Direct numerical studies of the correlation function of the dust support this point of view. 5. Discussion We have studied several models of surface growth in the quenched dis- ordered media near the pinning threshold. These models are in the universality class of the diode-resistor reverse percolation. Numerical results are in good agreement with anomalously large values of the critical exponents, obtained in many experiments in d = 1 + 1 and d = 2 + 1. The importance of quenched noise and pinning as a mechanism of surface roughening was suggested by several authors [26,27,31,32] which studied the continuum Langevin equations with quenched noise as models of surface growth. Our numerical results are in rather good agreement with the numerical results of those authors in d = 1 + 1. The differences in the values of exponents (less than 10%) may be explained by correction to scaling. Models A, B and C that we have studied here have strong corrections to scaling, but for large system sizes and times approach the same asymptotic behavior#’ . *I In d = 1 + 1, these correction vanishes for the systems of the order 214. In d = 2 + 1 and d = 3 + 1, they become negligible when L = 29 and L = 2’, respectively. In d = 4 + 1, they are present even for L = 25, which is the upper limit in our computational capacities. 210 S.V. Buldyrev et al. / Anisotropic percolation and surface roughening Acknowledgements We wish to thank T. Vicsek, A.L. Barabasi and G. Huber for helpful discussions, J. KertCsz for stimulating collaborations on some of the ideas presented here, D. Kessler for critical comments on the manuscript, and NSF for financial support. References [l] R. Jullien, J. Kerttsz, P. Meakin and D. Wolf, eds., Surface Disordering: Growth, Roughening, and Phase Transitions (Nova Science, New York, 1992). [2] T. Vicsek, Fractal Growth Phenomena, 2nd ed. (World Scientific, Singapore, 1992). [3] J. Krug and H. Spohn, in: Solids Far From Equilibrium: Growth, Morphology and Defects, C. Godreche, ed. (Cambridge Univ. Press, Cambridge, 1991). [4] P. Meakin, Phys. Rep., in press. [.5] T. Halpin-Healey and Y.-C. Zhang, Phys. Rep., in press. [6] M. Kardar, G. Parisi and Y.-C. Zhang, Phys. Rev. Lett. 56 (1986) 889; E. Medina, T. Hwa, M. Kardar and Y.-C. Zhang, Phys. Rev. A 39 (1989) 3053. [7] N. Martys, M. Cieplak and M.O. Robbins, Phys. Rev. Lett. 66 (1991) 1058. [8] J.M. Kim and J.M. Kosterlitz, Phys. Rev. Lett. 62 (1989) 2289. [9] B.M. Forrest and L.-H. Tang, Phys. Rev. Lett. 64 (1990) 1405. [lo] K. Moser, J. Kertesz and D.E. Wolf, Physica A 178 (1991) 215. [ll] T. Vicsek, M. Cserzo and V.K. Horvath, Physica A 167 (1990) 315. [12] M.A. Rubio, C.A. Edwards, A. Dougherty and J.P. Gollub, Phys. Rev. Lett. 63 (1990) 1685; V.K. Horvlth, F. Family and T. Vicsek, J. Phys. A 24 (1991) L25. [13] S.V. Buldyrev, A.-L. Barabasi, F. Caserta, S. Havlin, H.E. Stanley and T. Vicsek, Phys. Rev. A 45 (1992) R-8313. [14] S. Havlin, A.-L. Barabasi, S.V. Buldyrev, C.K. Peng, M. Schwartz, H.E. Stanley and T. Vicsek, in: Growth Patterns in Physical Sciences and Biology, E. Louis, L. Sander and P. Meakin, eds. (Plenum, New York, 1992). [15] A.-L. Barabasi, S.V. Buldyrev, S. Havlin, G. Huber, E. Stanley and T. Vicsek, in: Surface Disordering: Growth, Roughening, and Phase Transitions, R. Jullien, J. Kertesz, P. Meakin and D. Wolf, eds. (Nova Science, New York, 1992). [16] J. Zhang, Y.-C. Zhang, P. Alstrom and M.T. Levinsen, in: Surface Disordering: Growth, Roughening, and Phase Transitions, R. Jullien, J. Kertesz, P. Meakin and D. Wolf, eds. (Nova Science, New York, 1992). [17] M. Matsushita and S. Ouchi, Physica D 38 (1989) 246. [18] G. Dietler and Y.-C. Zhang, Physica A 191 (1992) 213. [19] S.V. Buldyrev, A.-L. Barabasi, S. Havlin, J. Kertesz, H.E. Stanley and H.S. Xenias, Physica A 191 (1992) 220. [20] J. Krim, I. Heyvaert, C. van Haesendonck and Y. Bruyhseraede, Phys. Rev. Lett. 70 (1993) 57. [21] L.-H. Tang and H. Leschhorn, Phys. Rev. A 45 (1992) R-8309. [22] D. Dhar, M. Barma, M.K. Phani, Phys. Rev. Lett. 47 (1981) 1238. [23] S. Redner, Phys. Rev. B 25 (1982) 3242. [24] S. Redner, in: Percolation Structures and Processes, G. Deutscher, R. Zallen and J. Adler, eds. (Hilger, Bristol, 1983). [25] S.V. Buldyrev, S. Havlin, J. Kerttsz, A. Shehter and H.E. Stanley, preprint. [26] G. Parisi, Europhys. Lett. 17 (1992) 673.  S.V. Buldyrev et al. 1 Anisotropic percolation and surface roughening 211 [27] Z. Csahbk, K. Honda, E. Somfai, M. Vicsek and T. Vicsek, Physica A 200 (1993) 136, these Proceedings. [28] W. Kinzel, in: Percolation Structures and Processes, G. Deutscher, R. Zallen and J. Adler, eds. (Hilger, Bristol, 1983). [29] A. Bunde and S. Havlin, eds., Fractals and Disordered Systems (Springer, Heidelberg, 1991). [30] J. KertCsz and H.J. Herrmann, J. Phys. A 18 (1985) L1109. [31] D.A. Kessler, H. Levine and Y. Tu, Phys. Rev. A 43 (1991) 4551. [32] M.H. Jensen and I. Procaccia, preprint.