Crossover from isotropic to directed percolation
We generalize the directed percolation (DP) model by relaxing the strict directionality of DP such that propagation can occur in either direction but with anisotropic probabilities. We denote the probabilities as and , with representing the average occupation probability and controlling the anisotropy. The Leath-Alexandrowicz method is used to grow a cluster from an active seed site. We call this model with two main growth directions biased directed percolation (BDP). Standard isotropic percolation (IP) and DP are the two limiting cases of the BDP model, corresponding to and respectively. In this work, besides IP and DP, we also consider the region. Extensive Monte Carlo simulations are carried out on the square and the simple-cubic lattices, and the numerical data are analyzed by finite-size scaling. We locate the percolation thresholds of the BDP model for and , and determine various critical exponents. These exponents are found to be consistent with those for standard DP. We also determine the renormalization exponent associated with the asymmetric perturbation due to near IP, and confirm that such an asymmetric scaling field is relevant at IP.
Directed percolation (DP), introduced in 1957 by Broadbent and Hammersley DP (), is a fundamental model in non-equilibrium statistical mechanics and represents the most common dynamic universality class Marro-Dickman (). DP has a very wide application, including flow in a porous rock in a gravitational field, forest fires, epidemic spreading, and surface chemical reactions Grassberger-JSP-79-13 (). The DP process can be illustrated in the simple example of bond DP on the square lattice. Along the horizontal (vertical) edges of the lattice, the propagation occurs in a particular direction only, e.g., toward the right (the up). Frequently, the preferred spreading direction is termed “temporal,” and the perpendicular one is called “spatial;” the two-dimensional DP is thus often called “(1+1)-dimensional DP.” The DP process has two distinct phases: the inactive phase for small occupation probability where the propagation quickly dies out, and the active phase for large . Between these two phases, a transition occurs at . As the threshold is approached, the temporal () and the spatial () correlation lengths diverge but with distinct critical exponents: and . The anisotropy is characterized by the so-called dynamic exponent . For , the order parameter , defined as the probability that a randomly selected site can generate an infinite cluster, becomes non-zero and its behavior can be described as , with another critical exponent. Below the upper critical dimensionality with , the three independent critical exponents, , , and , are sufficient to describe the DP universality class. While analytical results are scarce for DP, even in dimensions, approximation techniques like series expansion SE1 (); SE2 (); SE3 (); SE4 (); SE5 () and Monte Carlo simulations Grassberger-JPA-1989 (); Grassberger-Zhang (); Lubeck (); Voigt-Ziff () have produced fruitful results. Moreover, after a great deal of efforts, experimental realization of the DP process has been achieved Takeuchi (); errotum-Takeuchi (); extension-Takeuchi () in nematic liquid crystals, where the DP transition occurs between two turbulent states.
Analogously, standard isotropic percolation (IP) IP () is a fundamental model in equilibrium statistical mechanics. IP has attracted extensive research attention both in the physical and the mathematical communities, and the critical behavior is now well understood. Due to the isotropy, there exists only one spatial correlation length, which scales as near . Numerous exact results are now available in two dimensions (2D). For bond IP on the square lattice, the self-duality yields the threshold Kesten (); the values of are also exactly known for bond and site percolation on several other lattices Ziff-06-0 (); Ziff-06-1 (); Ziff-10-0 (), or have been determined to a high precision Feng08 (). Thanks to conformal field theory and Coulomb gas theory Cardy (); Smirnov (); Lawler (); Kesten-1987 (), the critical exponents and are also exactly known as and .
In this work, we introduce a generalized percolation propagation process that contains DP and IP as two special cases. On a given lattice, each edge is assigned to one of the three possible states: occupied by a directed bond along a particular direction, occupied by a directed bond against the particular direction, or unoccupied. This is illustrated in Fig. 1 (a), and the associated probabilities are denoted as , , and , respectively. As a result, the percolation process has two main growth directions. For , the symmetry between the two opposite directions is restored, and the system reduces to standard bond IP. In the limiting case or , propagation against or along the particular direction is forbidden, and one has standard DP. We call this percolation model with two main growth directions biased directed percolation (BDP). We note that the BDP model is described by the field-theoretic equation in Frey ().
A natural question arises: in between standard DP and IP, what is the nature of phase transition for BDP? For later convenience, we replace parameters and by two new variables as
The parameter is the average bond-occupation probability (irrespective of the bond direction), and accounts for the anisotropy. DP corresponds to or 1, while is for IP.
In this work, extensive Monte Carlo simulations are carried out for BDP in two and three dimensions. A dimensionless ratio is defined to locate the percolation threshold. The data are analyzed by finite-size scaling, and the critical exponents are determined. The numerical results suggest that the asymmetric perturbation due to is relevant near IP, and thus that as long as , BDP is in the DP universality class. These results further raise the following questions, remaining to be explored. For IP, is the asymmetric renormalization exponent a “new” critical exponent or related in some way to the known ones like and ? Particularly, can this “new” exponent be exactly obtained in two dimensions? If so, what is the exact value?
The remainder of this work is organized as follows. Section II introduces the BDP model, the sampled quantities, and the associated scaling behavior. Numerical results are presented in Secs. III and IV. A brief discussion is given in Sec. V.
Ii Model, sampled quantities, scaling behavior
We shall describe in details the BDP model on the square lattice. The generalization to higher dimensions is straightforward.
As usual in the study of DP or IP, we view the BDP model as a stochastic growth process, and use the Leath-Alexandrowicz method Leath (); Alex () to grow the percolation cluster starting from a seed site. Given the square lattice and the seed “o” in Fig. 1(b), for each of the neighboring edges of site o, a random number is drawn to determine the edge state. If and only if the edge is occupied and the direction originates from the seed o, the neighboring site is activated and belongs to the growing cluster. For instance, in Fig. 1(b), the four neighboring edges of site o are all occupied, but site c remains inactivated because of the “wrong” direction. After all the four neighboring edges have been visited, one continues the growing procedure from the newly added sites. In other words, one grows the percolation cluster shell by shell (the breadth-first scheme). The growth of the cluster continues until the procedure dies out or the maximum distance is reached, which is set at the beginning of the simulation.
ii.2 Sampled quantities
In the cluster-growing process, the number of activated sites is recorded as a function of the shell number . Let us count the shell of site “o” to be the first shell, the configuration in Fig. 1(b) has for , respectively. Besides , one also records the Euclidean distance of each activated site to the seed “o” for IP and to the axis for the anisotropic case. The reason for using different definitions of is that, for the anisotropic case, the average center of activated sites is expected to move linearly along the preferred direction, as increases. Accordingly, we define a revised gyration radius as
The statistical averages and are then measured, as well as their statistical uncertainties. We also measure the survival probability that at least one site remains activated at the th shell and the accumulated activated site number .
In Monte Carlo study of critical phenomena and phase transitions, it is found that dimensionless ratios like the Binder cumulant are very useful in locating the critical point. Therefore, we also define a dimensionless ratio .
ii.3 Scaling behavior
Near the percolation threshold , one expects the following scaling behavior
where represents a small deviation from . Symbols , and denote the associated critical exponents, and , , and are universal functions. For simplicity, only one scaling field, which accounts for the effect due to deviation from , is explicitly included in Eq. (3). Right at , as increases, the survival probability decays to zero while the other quantities diverge, except for the ratio which goes to a constant. A trivial relation is .
For standard DP ( or 1), exponents and are normally denoted as and , respectively ( is also denoted as in Hinrichsen-2000 ()). It can be shown that exponent is . Further, exponent relates to and the dynamic exponent as , where arises from the behavior . Below the upper critical dimensionality with , there exist three independent exponents, which can be chosen as , , and . The others can be obtained by the scaling relations Hinrichsen-2000 ()
where the last one involves the spatial dimensionality and is called the hyperscaling relation. In dimensions, these exponents have been determined to high precision: , , and SE4 (). In dimensions, these exponents are , , and Grassberger-Zhang (); Voigt-Ziff (); Perlsman (); Junfeng ().
For standard IP (), the shell number is frequently called “chemical distance” Havlin (), accounting for the minimum length among all the possible paths between the seed site and the activated sites on the th shell. At , the length of the chemical path relates to the Euclidean distance as Grassberger-JPA-1992 (); Grassberger-JPA-1985 (), with denoting the shortest-path exponent. In terms of the Euclidean distance , it is known that the survival probability scales as , the accumulated site number , and the -deviating scaling behavior . This immediately yields , and . For IP, one has the scaling relation as
In 2D, and are exactly known as and , which yield and . The shortest-path exponent , together with the so-called backbone exponent, is among the few critical exponents of which the exact values are not known for the 2D percolation universality class. It was conjectured to be Deng-PRE-81 (), and some recent estimates are 1.1306(3) Grassberger-JPA-32 () and 1.13078(5) JY (). In three dimensions, no exact results are available, and the numerical estimates are Grassberger-JPA-25 (), and deng (), which yield .
In this work, we consider the BDP model on the square lattice for 2D and the simple-cubic lattice for 3D. The simulation applies the aforementioned Leath-Alexandrowicz growth method. The dimensionless ratio is used to locate the percolation threshold . According to Eq. (3), ratio is expected to have an approximate common intersection at for different shell number . At the threshold , as , the common intersection converges to a universal value and the slope of increases as .
iii.1 Standard IP
Standard IP corresponds to . Monte Carlo simulation was carried out up to for 2D and for 3D. About samples were taken for each data point on each lattice. The data are shown in Fig. 2. Indeed, we find an approximate common intersection near and for 2D and 3D, respectively. This agrees with the known threshold and Lorenz (). Note that, since the occupied bond can propagate the growth process only if it has the correct orientation, there is a factor of 2 difference between the bond-occupation probability here and the of the equivalent bond percolation probability.
To have a better estimate of , according to a least-squared criterion, the data are fitted by
which is obtained by Taylor-expanding Eq. (3) and taking into account finite-size corrections due to the leading irrelevant scaling field and analytical background contribution. These are described by the two terms with amplitudes and , of which the term with arises from the nonlinearity of the relevant scaling field in terms of the deviation , and the one with accounts for the combined effect of the leading relevant and irrelevant scaling fields. In the fits, various formulas are tried, which correspond to different combinations of those terms in Eq. (6). For a given formula, the data for small are gradually excluded from the fits to see how the residual changes with respect to . The results from different formulas are compared with each other to estimate the possible systematic errors. In two dimensions, we obtain , , , and . Note that the leading irrelevant thermal scaling field is for 2D percolation universality Ziff-PRE-2011 (); apparently, the leading correction exponent does not correspond to . Instead, should be associated with the chemical distance. From the relations , , and , and the exact values and , we determine from , and from .
In three dimensions, our results are , , , and . Our estimate of agrees with the existing one Lorenz (), and has a comparable error margin.
|(known)||IP (); Cardy (); Smirnov (); Lawler (); Kesten-1987 ()||IP (); Cardy (); Smirnov (); Lawler (); Kesten-1987 ()||Grassberger-JPA-32 (); Deng-PRE-81 (); JY ()||IP (); Kesten ()|
|(known)||deng ()||Grassberger-JPA-25 ()||Lorenz ()|
To estimate other critical exponents, we simulate right at the threshold for 2D and for 3D. The simulation was carried out for up to for 2D and 2048 for 3D. Further, to eliminate one more unknown parameter in the fits, we measure the dimensionless ratios and . These data are fitted by
In two dimensions, the results are and for , and and for . For all these three ratios, the leading correction is described by an exponent . Taking into account the exact values , one has from and from .
In three dimensions, the results are and for , and and for . Combining the estimate and together, one has and . Our result agrees well with the existing result Grassberger-JPA-25 (), and significantly improves the error margin.
For comparison, these results are summarized in Table 1.
iii.2 Standard DP
We simulate standard DP by taking . The simulation was carried out for up to for 2D, and 2048 for 3D. The number of samples for each data point is about in 2D and in 3D.
The data are shown in Fig. 3. A good intersection is observed for both 2D and 3D, which yields for 2D and for 3D, from a rough visual fitting.
We fit the data more precisely using Eq. (6). On the square lattice, we obtain , , , and . The estimate of the percolation threshold agrees well with the existing more precise result SE4 (). From the relations and , we have and . On the simple-cubic lattice, our results are , , , which yield and . Here the is too small to estimate since the numerical data of can be well described even though we do not include any corrections. The agreement of with the existing estimate Junfeng () is within two standard deviations.
Analogously, we simulate right at the percolation threshold for 2D and for 3D. The dimensionless ratios and are measured, and the data are fitted by Eq. (7). For 2D, the results are , and , , which yield and . Taking into account the estimates of and , one has . For 3D, the results are and , which yield that and .
These results are listed in Table 2.
For the purpose of studying BDP, we choose and . The simulation was carried out for up to for 2D and 2048 for 3D. About samples were taken for each data point in each case.
The data are shown in Fig. 4 for 2D and Fig. 5 for 3D. The transitions are also clearly observed, but the approximate common intersections are not as good as those for standard DP and IP. This suggests the existence of additional finite-size corrections.
The data are also fitted by Eq. (6) according to a least-squared criterion. To account for the possible existence of additional corrections, we replace the terms in Eq. (6), with , , and , by . The exponent is fixed at , in accordance with our above estimate of for both standard IP and DP. Indeed, the new source of finite-size correction can be identified in the fits, which yield both in 2D and 3D. The results for , , and are summarized in Table 2.
|1||0.644 700 5(8)||0.277(2)||1.736(9)||1.580 6(3)||0.313 7(2)||0.159 44(9)|
|0.8||0.768 708(1)||0.278(2)||1.74(1)||1.577(5)||0.314 1(4)||0.159 5(1)|
|0.6||0.929 668(3)||0.279(2)||1.754(6)||1.578(5)||0.316 1(8)||0.159(1)|
|3||Junfeng ()||1||0.382 224 64(4)||0.581 2(6)||1.289 0(7)||1.766 5(2)||0.230 81(7)||0.450 9(2)|
|1||0.382 225 6(5)||0.582(5)||1.287(4)||1.767(3)||0.231 2(1)||0.451 9(8)|
The determination of the critical exponents and is obtained in an analogous way by simulating at the estimated percolation threshold, and the results are listed in Table 2.
The results in Table 2 strongly suggest that, as long as deviates from , the system falls into the standard DP universality class. For an illustration, we make the log-log plot of the critical quantity versus the shell number in Fig. 6. Clearly, the slope for is distinct from those for the other cases, which are independent of .
Iv Crossover exponent
The fact that BDP for is in the DP universality means that in the language of renormalization group theory, the operator associated with the asymmetric perturbation is relevant near the IP fixed point. To confirm this, we simulate BDP near IP with by varying . The simulation is up to , and is set at , and . The results for in two dimensions are shown in Fig. 7 versus ; note that BDPs for are identical. These data are also analyzed by Eq. (6) with being replaced by the exponent for the symmetric scaling field and the odd terms with respect to being set zero. We obtain , which suggests that may be exactly .
According to scaling theory, the phase transition line approaches to the critical IP as Riedel-1969 ()
where is the so-called crossover exponent. We carried out some Monte Carlo simulations and determined a set of critical points near IP; they are 13 critical points with . In Fig. 7, we plot versus in log-log scale, which indeed has slope approximately equal to .
We also perform a similar study near the critical IP in 3D, and obtain and .
We introduce a biased directed percolation model, which includes standard isotropic and directed percolation as two special cases. Large-scale Monte Carlo simulations are carried out in two and three dimensions. We find that the operator associated with the anisotropy is relevant near the IP fixed point, which implies that BDP in the region is in the DP university class. On this basis, the phase diagram and the associated renormalization flows are shown in Fig. 8.
Since the upper critical dimensionality is different for standard IP and DP, it is not clear whether the similar renormalization flows would hold in higher dimensions. We mention that such crossover phenomena have attract much attention both in the fields of equilibrium and non-equilibrium statistical mechanics Pfeuty (); Aharony (); Lubeck-JSM (); Mendes (); Frojdh (); Janssen (); Schonmann-JSP-1986 (). In retrospect, it is not surprising that the asymmetric perturbation is relevant near IP. At IP, all the directions are equivalent and “spatial” and “temporal” directions cannot be defined. However, as soon as , such a symmetry is broken and the center of the activated sites moves along the “temporal” direction as the growing process continues. It is also plausible that as long as the “spatial” and “temporal” symmetry is not restored, such an asymmetric perturbation is irrelevant near DP. This is similar to the fact that asymmetric diffusion on the basic contact process is irrelevant Schonmann-JSP-1986 (). In terms of the chemical distance , the effect from the anisotropy can be asymptotically described as with and . One can also use the Euclidean distance to describe such an anisotropic effect as with . Substituting the value into , one obtains and .
When viewing standard isotropic percolation in the framework of BDP, one observes that two independent critical exponents, e.g., and , are no longer sufficient to describe the critical scaling behavior. In this case, the shortest-path exponent appears naturally and becomes indispensable, and thus isotropic percolation also has three independent critical exponents. Our estimate of significantly improves over the existing results both in two and three dimensions. Our result does not agree with the recently conjectured value Deng-PRE-81 () in two dimensions. This result appears to refute the conjectured value. On the other hand, we note that, in terms of the chemical distance , a new source of finite-size corrections occurs in the scaling behavior, and these corrections are not well understood yet. Further, we observe that the restored symmetry for IP can be regarded as in the BDP model. In some cases, the coincidence of two critical exponents may suggest the existence of logarithmic corrections of the or form, and they can be either additive or multiplicative. In practice, logarithmic finite-size corrections have indeed been observed for standard isotropic percolation in two dimensions Feng08 (), which is in terms of Euclidean distance. In this sense, we cannot entirely exclude the possibility that the tiny difference between the present numerical result for and the conjectured value arises from some unknown corrections that have not been taken into account in the numerical analysis. Numerical investigation of this problem seems very difficult if not impossible. Nevertheless, since the exact value of is conjectured as a function of for the -state Potts model Deng-PRE-81 (), one can accumulate more numerical evidence by studying the case.
Finally, the numerical estimate of the critical exponent or due to the asymmetric perturbation near IP raises a question: is it a “new” independent critical exponent or simply related in some way to the known ones like , , and ? In particular, in two dimensions, one would ask whether or can be exactly obtained in the framework of Stochastic Loewner Evolution (SLE), conformal field theory or Coulomb gas theory.
This work was supported in part by NSFC under Grant No. 10975127 and 91024026, and the Chinese Academy of Science. RMZ acknowledges support from National Science Foundation Grant No. DMS-0553487. We also would like to thank Dr. Timothy M. Garoni in Monash University for valuable comments.
- (1) S. R. Broadbent and J. M. Hammersley, Proc. Camb. Phil. Soc. 53, 629 (1957).
- (2) J. Marro and R. Dickman, Nonequilibrium phase transitions in lattice models, Cambridge University Press, Cambridge, 1999.
- (3) P. Grassberger, J. Stat. Phys. 79, 13 (1995).
- (4) K. De’Bell, and J. W. Essam, J. Phys. A 16, 3553 (1983).
- (5) J. W. Essam, K. De’Bell, J. Adler, F. M. Bhatti, Phys. Rev. B 33, 1982 (1986).
- (6) R. J. Baxter and A. J. Guttmann, J. Phys. A 21, 3193 (1988).
- (7) I. Jensen, J. Phys. A 32, 5233 (1999).
- (8) I. Jensen, J. Phys. A 37, 6899 (2004).
- (9) P. Grassberger, J. Phys. A 22, 3673 (1989).
- (10) P. Grassberger and Y. C. Zhang, Physica A 224, 169 (1996).
- (11) C. A. Voigt and R. M. Ziff, Phy. Rev. E 56, R6241 (1997).
- (12) S. Lübeck and R. D. Willmam, J. Stat. Phys, 115, 516 (2004).
- (13) K. A. Takeuchi, M. Kuroda, H. Chatè and M. Sano, Phys. Rev. Lett. 99, 234503 (2007).
- (14) K. A. Takeuchi, M. Kuroda, H. Chatè and M. Sano, Phys. Rev. Lett. 103, 089901 (E) (2009).
- (15) K. A. Takeuchi, M. Kuroda, H. Chatè and M. Sano, Phys. Rev. E 80, 051116 (2009).
- (16) D. Stauffer and A. Aharony, Introduction to Percolation Theory (Taylor and Francis, London, 1992), and references therein.
- (17) H. Kesten, Comm. Math. Phys. 74, 41 (1980).
- (18) R. M. Ziff, Phys. Rev. E 73, 016134 (2006).
- (19) R. M. Ziff and C. R. Scullard, J. Phys. A 39, 15083 (2006).
- (20) C. R. Scullard and R. M. Ziff, J. Stat. Mech. P03021 (2010).
- (21) X. Feng, Y. Deng and H. W. J. Blöte, Phys. Rev. E 78, 031136 (2008).
- (22) J. L. Cardy, Nucl. Phys. B 240, 514 (1984).
- (23) S. Smirnov, W. Werner, Math. Research Letters, 8, 729 (2001).
- (24) G. F. Lawler, O. Schramm, W. Werner, One-arm exponent for critical 2D percolation. Electron. J. Probab. 7, paper no.2 (2002).
- (25) H. Kesten, Comm. Math. Phys. 109, 109 (1987).
- (26) E. Frey, U. C. Taüber, F. Schwabl, Phys. Rev. E 49, 5058 (1994).
- (27) P. L. Leath, Phys. Rev. B 14, 5046 (1976).
- (28) Z. Alexandrowicz, Phys. Lett. A 80, 284 (1980).
- (29) H. Hinrichsen, Adv. Phys. 49, 815 (2000).
- (30) E. Perlsman and S. Havlin, Europhys. Lett. 58, 176 (2002).
- (31) J. Wang, Q. Liu and Y. Deng, arXiv:1201.3006.
- (32) S. Havlin, B. Trus, G. H. Weiss and D. ben-Avraham, J. Phys. A 18, L247 (1985).
- (33) P. Grassberger, J. Phys. A 25, 5475 (1992).
- (34) P. Grassberger, J. Phys. A 18, L215 (1985).
- (35) Y. Deng, W. Zhang, T. M. Garoni, A. D. Sokal, and A. Sportiello, Phys. Rev. E 81, 020102(R) (2010).
- (36) P. Grassberger, J. Phys. A 32, 6233 (1999).
- (37) J. Yang, Z. Zhou, R. M. Ziff and Y. Deng, in preparation (2012).
- (38) P. Grassberger, J. Phys. A 25, 5867 (1992).
- (39) Y. Deng and H. W. J. Blöte, Phys. Rev. E 72, 016126 (2005).
- (40) C. D. Lorenz and R. M. Ziff, Phys. Rev. E 57, 230 (1998).
- (41) R. M. Ziff, Phys. Rev. E 83, 020107 (2011).
- (42) E. K. Riedel and F. J. Wegner, Z. Phys. 225, 195 (1969).
- (43) P. Pfeuty, G. Toulouse, Introduction to the renormalization group and critical phenomena (John Wiley & Sons, Chichester,1994).
- (44) A. Aharony, Dependence of universal critical behavior on symmetry and range of interaction in Phase Transition and Critical Phenomena, Vol.6, edited by C. Domb and M.S. Green (Academic Press, London, 1976).
- (45) S. Lübeck, J. Stat. Mech. P09009 (2006).
- (46) J. F. F. Mendes, R. Dickman, H. Herrmann, Phys. Rev. E 54, R3071 (1996).
- (47) P. Fröjdh and M. den Nijs, Phys. Rev. Lett. 78, 1850 (1997).
- (48) H. K. Janssen, O. Stenull, Phys. Rev. E 62, 3173 (2000).
- (49) R. H. Schonmann, J. Stat. Phys. 44, 505 (1986).