- Email: [email protected]
Universal scaling functions of a dissipative Manna model
Computer Physics Communications 182 (2011) 226–228
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Universal scaling functions of a dissipative Manna model Chien-Fu Chen, An-Chung Cheng, Yi-Duen Wang, Chai-Yu Lin ∗ Department of Physics, National Chung Cheng University, Chia-Yi 621, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 1 March 2010 Received in revised form 19 June 2010 Accepted 16 July 2010 Available online 29 July 2010
This study uses a losing probability f to measure the magnitude of the bulk dissipation of a dissipative Manna model where P (s; r , L , f ), the probability distribution of toppling size s, on a L × ( L /r ) square lattice is calculated. This approach assumes that the bulk dissipation corresponds to a characteristic length ξ . The finding ξ ∼ f −ν leads to P (s; r , L , f ) = s−τ H (s/[ L / A r ] D , [ L / A r ] f ν ) where H is a twovariable universal scaling function, τ , D, and ν are exponents, and A r is a nonuniversal metric factor determined by the calculation of moment. © 2010 Elsevier B.V. All rights reserved.
Keywords: Manna model Universal scaling function Bulk dissipation
Sandpile models have been proposed to explain the concept of self-organized criticality (SOC) [1,2]. The balance of external driving flow and dissipation flow maintains balance in a sandpile system. Several researchers [3–5] have pointed out that the criticality of a sandpile model occurs when (i) the driving rate h is zero, (ii) the dissipation rate is zero, and (iii) the ratio of h to is zero. Some specific sandpile models adopt absorbing phase transition models [6] and linear interface depinning models [7] to increase the precise determination of universality class. In addition, the renormalization group approach [8–10] leads to a much better understanding of the insight of the sandpile models. Moment analysis is a powerful tool for numerically analyzing the exponents of a specific sandpile model. Using the moment analysis, researchers have found that the Manna sandpile [11] exhibits finite-size scaling behavior [12,13], while the Bak–Tang–Wiesenfeld sandpile [1] exhibits multiscaling behavior [14]. A sandpile model is built on a lattice with the open boundary condition or periodic boundary condition. In this case, dissipation occurs through the boundary or the bulk. The boundary (bulk) dissipation rate L (b ) describes the magnitude of the boundary (bulk) dissipation where = L + b . The periodic boundary condition leads to L = 0. For the open boundary condition, the value of L depends on the lattice size L where L → 0 as L → ∞. There is an interesting connection [15,16] between the Abelian variant of the Manna sandpile models and random-walk models. Using a L × L square lattice, previous researchers have accurately determined L = (0.07L + 0.28L 2 )−1 . The terms L and b control the avalanche size s of a sandpile model [17]. Previous research [16] shows that the behavior
*
Corresponding author. E-mail address: [email protected] (C.-Y. Lin).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.07.033
© 2010
Elsevier B.V. All rights reserved.
of a dissipative Manna model with ( L , b ) = (a, 0) is similar to one with ( L , b ) = (0, a) where a is a constant. Thus, b and L play similar role. Since L corresponds to a length L, b should correspond to the length ξ , called the characteristic length. This paper considers a dissipative variant of the Abelian Manna twostate model in which L = 0 and b = 0. Compared with some original SOC models, which only consider boundary dissipation, this model reveals the competition between L and ξ and provides an insight into the scaling behaviors. This study defines a dissipative variant of the Abelian Manna two-state model [11–13] on a L × ( L /r ) square lattice with open boundary condition, where r is the aspect ratio. In this case, each lattice site is associated with a nonnegative integer as its height. A site topples when its height is equal to or greater than zc where zc = 2 is the toppling threshold. The toppling rule is to (i) decrease the height of the toppling site by zc , and (ii) operate the transferring procedure zc times. The transferring procedure increases the height of one randomly chosen nearest neighbor of the toppling site by one (or zero) with a probability of 1 − f (or f ). This sandpile system exhibits two methods of dissipation: (1) the mechanism of vanishing energy, called bulk dissipation, where b = zc f , and (2) the mechanism of transferring energy to out of lattice, called boundary dissipation. After driving the system by adding one unit of energy to a randomly chosen site, the number of topplings s determines the size of the avalanche. Repeating the adding and relaxing procedures many times generates many avalanches. Every simulation in this paper conducts the adding procedure at least 2.2 × 107 times, and then excludes the data of the first 2 × 106 times from statistical analysis. According to Ref. [16], a characteristic length ξ describes the behavior of the system. This study sets ξ ∼ b−ν ∼ f −ν , where ν = 0.5 is the characteristic length exponent [17]. Consider P (s; r , L , f ), which is the probability distribution of s for a fixed f
C.-F. Chen et al. / Computer Physics Communications 182 (2011) 226–228
Fig. 2. sτ P (s; r , L , f ) as a function of s/ L D .
Fig. 1. P (s; r , L , f ) as a function of s.
and a fixed L. Based on the general form of a scaling function [18], P (s; r , L , f ) can be expressed as
P (s; r , L , f ) = s−τ G s/ L D , L /ξ, r ,
(1)
where G depending on r is a three-variable scaling function and the condition ξ L (L ξ ) represents to critical (finite-size) scaling. According to L /ξ ∼ L f ν , this study simulates P (s; r , L , f ) for sets ir , iir , and iiir where r denotes the aspect ratio of one specific set. For these sets, r = 1, 2, and 3. These sets correspond to L f ν = 109.39 (i3 ), 79.21 (i2 ), 50 (i1 ), 21.88 (ii3 ), 15.85 (ii2 ), 10 (ii1 ), 3.46 (iii3 ), 2.50 (iii2 ), and 1.58 (iii1 ). For each set, L = 600 and 300 are chosen. Thus, the simulations produce 18 values of (r , L , f ). The simulation parameters are ( L , f ) = (600, 0.03324) and (300, 0.133162) for set i3 , (600, 0.01743) and (300, 0.0698336) for set i2 , (600, 0.00694) and (300, 0.0278) for set i1 , (600, 0.00133) and (300, 0.005269) for set ii3 , (600, 0.000698) and (300, 0.0027632) for set ii2 , (600, 0.000278) and (300, 0.0011) for set ii1 , (600, 0.0000332) and (300, 0.000132683) for set iii3 , (600, 0.00001741) and (300, 0.0000695824) for set iii2 , and (600, 0.00000693) and (300, 0.0000277) for set iii1 . Figs. 1 and 2 plot P (s; r , L , f ) and sτ P (s; r , L , f ) as a function of s and s/ L D , respectively, where τ = 1.27 and D = 2.75 [13]. The 18 curves in Fig. 1 almost fall on the nine collapsed curves shown in Fig. 2 where each coalesced curve corresponds to a value of L f ν and a fixed r. These simulations suggest the existence of the scaling function G. Assume that the characteristic length of a system with aspect ratio r satisfies
ξ(r , f ) = A r f −ν ,
227
(2)
where A r is a nonuniversal metric factor for a given r. Consider s = sP (s; r , L , f )ds and simulate two cases: (a) f = 0 and a given finite L (i.e., L ξ ), which corresponds to boundary dissipation dominating the system. In this case, L determines system behavior. (b) L → ∞ and a given finite f > 0 (i.e., ξ L), which corresponds to bulk dissipation dominating the system. In this case, ξ determines system behavior. For case (a), consider r = 1, 2, and 3. Simulating P (s; r , L , f = 0) for L = 300, 600, 900, and 1200 shows that s = cr L σ . For case (b), choose a large L to represent the condition L → ∞. Then, select L = 1200 and set r = 1, 2, and 3. Simulating P (s; r , L , f ) for f = 10−1 , 10−2 , 10−3 , and 10−4 shows that s = er f −β . Taken together, cases (a) and (b) imply that
s =
cr L σ
dr ξ σ = er f −β
for L ξ,
(3)
for ξ L .
Previous research shows σ = 1.97 [13]. The knowledge of critical phenomena suggests that σ = σ . Then, β = σ ν . Numerical simulations give that the terms cr , dr , and er are constants for a fixed r while er = e is independent of r. Assuming that σ σ
cr
c1
=
dr
d1
and
dr
d1
=
er /[ A r ]
e 1 /[ A 1 ]σ
=
A1 Ar
(4)
leads to the following relationship
Ar =
c1 cr
1/σ A1.
(5)
Since the value of A r / A 1 is important and the real vale of A 1 is irrelevant in this study, set A 1 = 1 for simplification leading to A r = [c 1 /cr ]1/σ . Thus, this study employs the values of cr obtained from the first moment analysis (i.e., s ) and Eq. (5) to calculate the values of A r . This leads to values of (cr , A r ) = (0.114, 1.000) for r = 1, (0.046, 1.582) for r = 2, and (0.024, 2.199) for r = 3 where the error bars are within 0.005. To unify G with different r into a universal scaling function, this study proposes
P (s; r , L , f ) = s−τ H s/[ L / A r ] D , [ L / A r ]/ f −ν ,
(6)
where H is a two-variable function. To verify Eq. (6), Fig. 3 plots sτ P as a function of s/( L / A r ) D . The 18 curves in Fig. 1 almost fall on the 3 coalesced curves shown in Fig. 3 where sets ir , iir , and iiir for r = 1, 2, and 3 fall on sets I, II, and III, respectively. Note that sets I, II, and III correspond to L /ξ = [ L / A r ]/ f −ν ) = 50, 10, and 1.58, respectively. This suggests the existence of the universal scaling function H . This study uses the critical phenomena framework to investigate a SOC model with different aspect ratios, where a losing probability f functions as a tuning parameter and f = 0 simulates the critical point of the system. Based on the connection between f and ξ , this study proposes the universal scaling function form. The proposed procedure is applicable to calculating the universal scaling functions of other physical quantities of sandpiles, e.g., the toppling area or duration of an avalanche. Finally, fitting exponents τ and D requires the preparations of several (L, f , A r ) which satisfy the condition that [ L / A r ] f ν is a constant according to the expression of Eq. (6). Otherwise, the obtained exponent values will be inaccurate.
228
C.-F. Chen et al. / Computer Physics Communications 182 (2011) 226–228
References
Fig. 3. sτ P (s; r , L , f ) as a function of s/[ L / A r ] D .
Acknowledgement C.-Y. Lin acknowledges the support of the National Science Council of Taiwan, R.O.C. under Grant No. NSC 97-2112-M-194003-MY2.
[1] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381; P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. A 38 (1988) 364. [2] H.J. Jensen, Self-Organized Criticality, Cambridge University Press, New York, 1998. [3] D. Sornette, A. Johansen, I. Dornic, J. Phys. I 5 (1995) 325. [4] A. Vespignani, S. Zapperi, V. Loreto, Phys. Rev. Lett. 77 (1996) 4560. [5] A. Vespignani, S. Zapperi, Phys. Rev. Lett. 78 (1997) 4793; A. Vespignani, S. Zapperi, Phys. Rev. E 57 (1998) 6345. [6] R. Dickman, A. Vespignani, S. Zapperi, Phys. Rev. E 57 (1998) 5095; A. Vespignani, R. Dickman, M.A. Muñoz, S. Zapperi, Phys. Rev. Lett. 81 (1998) 5676. [7] M. Paczuski, S. Boettcher, Phys. Rev. Lett. 77 (1996) 111. [8] L. Pietronero, A. Vespignani, S. Zapperi, Phys. Rev. Lett. 72 (1994) 1690. [9] E.V. Ivashkevich, Phys. Rev. Lett. 76 (1996) 3368. [10] C.-Y. Lin, A.-C. Cheng, T.-M. Liaw, Phys. Rev. E 76 (2007) 041114; A.-C. Cheng, C.-F. Chen, C.-Y. Lin, Int. J. Mod. Phys. C 19 (2008) 1695; C.-Y. Lin, Phys. Rev. E 81 (2010) 021112. [11] S.S. Manna, J. Phys. A 24 (1991) L363. [12] S. Lübeck, Phys. Rev. E 61 (2000) 204; S. Lübeck, Int. J. Mod. Phys. B 18 (2004) 3977. [13] A. Chessa, A. Vespignani, S. Zapperi, Comput. Phys. Commun. 121–122 (1999) 299. [14] C. Tebaldi, M. De Menech, A.L. Stella, Phys. Rev. Lett. 83 (1999) 3952. [15] Y. Shilo, O. Biham, Phys. Rev. E 67 (2003) 066102. [16] O. Malcai, Y. Shilo, O. Biham, Phys. Rev. E 73 (2006) 056125. [17] C.-Y. Lin, C.-F. Chen, C.-N. Chen, C.-S. Yang, I.-Min Jiang, Phys. Rev. E 74 (2006) 031304; C.-F. Chen, C.-Y. Lin, Int. J. Mod. Phys. C 20 (2009) 273. [18] H.E. Stanley, Introduction to Phase Transitions and Critical Phenomena, Oxford Univ. Press, New York, 1971.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Universal scaling functions of a dissipative Manna model Chien-Fu Chen, An-Chung Cheng, Yi-Duen Wang, Chai-Yu Lin ∗ Department of Physics, National Chung Cheng University, Chia-Yi 621, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 1 March 2010 Received in revised form 19 June 2010 Accepted 16 July 2010 Available online 29 July 2010
This study uses a losing probability f to measure the magnitude of the bulk dissipation of a dissipative Manna model where P (s; r , L , f ), the probability distribution of toppling size s, on a L × ( L /r ) square lattice is calculated. This approach assumes that the bulk dissipation corresponds to a characteristic length ξ . The finding ξ ∼ f −ν leads to P (s; r , L , f ) = s−τ H (s/[ L / A r ] D , [ L / A r ] f ν ) where H is a twovariable universal scaling function, τ , D, and ν are exponents, and A r is a nonuniversal metric factor determined by the calculation of moment. © 2010 Elsevier B.V. All rights reserved.
Keywords: Manna model Universal scaling function Bulk dissipation
Sandpile models have been proposed to explain the concept of self-organized criticality (SOC) [1,2]. The balance of external driving flow and dissipation flow maintains balance in a sandpile system. Several researchers [3–5] have pointed out that the criticality of a sandpile model occurs when (i) the driving rate h is zero, (ii) the dissipation rate is zero, and (iii) the ratio of h to is zero. Some specific sandpile models adopt absorbing phase transition models [6] and linear interface depinning models [7] to increase the precise determination of universality class. In addition, the renormalization group approach [8–10] leads to a much better understanding of the insight of the sandpile models. Moment analysis is a powerful tool for numerically analyzing the exponents of a specific sandpile model. Using the moment analysis, researchers have found that the Manna sandpile [11] exhibits finite-size scaling behavior [12,13], while the Bak–Tang–Wiesenfeld sandpile [1] exhibits multiscaling behavior [14]. A sandpile model is built on a lattice with the open boundary condition or periodic boundary condition. In this case, dissipation occurs through the boundary or the bulk. The boundary (bulk) dissipation rate L (b ) describes the magnitude of the boundary (bulk) dissipation where = L + b . The periodic boundary condition leads to L = 0. For the open boundary condition, the value of L depends on the lattice size L where L → 0 as L → ∞. There is an interesting connection [15,16] between the Abelian variant of the Manna sandpile models and random-walk models. Using a L × L square lattice, previous researchers have accurately determined L = (0.07L + 0.28L 2 )−1 . The terms L and b control the avalanche size s of a sandpile model [17]. Previous research [16] shows that the behavior
*
Corresponding author. E-mail address: [email protected] (C.-Y. Lin).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.07.033
© 2010
Elsevier B.V. All rights reserved.
of a dissipative Manna model with ( L , b ) = (a, 0) is similar to one with ( L , b ) = (0, a) where a is a constant. Thus, b and L play similar role. Since L corresponds to a length L, b should correspond to the length ξ , called the characteristic length. This paper considers a dissipative variant of the Abelian Manna twostate model in which L = 0 and b = 0. Compared with some original SOC models, which only consider boundary dissipation, this model reveals the competition between L and ξ and provides an insight into the scaling behaviors. This study defines a dissipative variant of the Abelian Manna two-state model [11–13] on a L × ( L /r ) square lattice with open boundary condition, where r is the aspect ratio. In this case, each lattice site is associated with a nonnegative integer as its height. A site topples when its height is equal to or greater than zc where zc = 2 is the toppling threshold. The toppling rule is to (i) decrease the height of the toppling site by zc , and (ii) operate the transferring procedure zc times. The transferring procedure increases the height of one randomly chosen nearest neighbor of the toppling site by one (or zero) with a probability of 1 − f (or f ). This sandpile system exhibits two methods of dissipation: (1) the mechanism of vanishing energy, called bulk dissipation, where b = zc f , and (2) the mechanism of transferring energy to out of lattice, called boundary dissipation. After driving the system by adding one unit of energy to a randomly chosen site, the number of topplings s determines the size of the avalanche. Repeating the adding and relaxing procedures many times generates many avalanches. Every simulation in this paper conducts the adding procedure at least 2.2 × 107 times, and then excludes the data of the first 2 × 106 times from statistical analysis. According to Ref. [16], a characteristic length ξ describes the behavior of the system. This study sets ξ ∼ b−ν ∼ f −ν , where ν = 0.5 is the characteristic length exponent [17]. Consider P (s; r , L , f ), which is the probability distribution of s for a fixed f
C.-F. Chen et al. / Computer Physics Communications 182 (2011) 226–228
Fig. 2. sτ P (s; r , L , f ) as a function of s/ L D .
Fig. 1. P (s; r , L , f ) as a function of s.
and a fixed L. Based on the general form of a scaling function [18], P (s; r , L , f ) can be expressed as
P (s; r , L , f ) = s−τ G s/ L D , L /ξ, r ,
(1)
where G depending on r is a three-variable scaling function and the condition ξ L (L ξ ) represents to critical (finite-size) scaling. According to L /ξ ∼ L f ν , this study simulates P (s; r , L , f ) for sets ir , iir , and iiir where r denotes the aspect ratio of one specific set. For these sets, r = 1, 2, and 3. These sets correspond to L f ν = 109.39 (i3 ), 79.21 (i2 ), 50 (i1 ), 21.88 (ii3 ), 15.85 (ii2 ), 10 (ii1 ), 3.46 (iii3 ), 2.50 (iii2 ), and 1.58 (iii1 ). For each set, L = 600 and 300 are chosen. Thus, the simulations produce 18 values of (r , L , f ). The simulation parameters are ( L , f ) = (600, 0.03324) and (300, 0.133162) for set i3 , (600, 0.01743) and (300, 0.0698336) for set i2 , (600, 0.00694) and (300, 0.0278) for set i1 , (600, 0.00133) and (300, 0.005269) for set ii3 , (600, 0.000698) and (300, 0.0027632) for set ii2 , (600, 0.000278) and (300, 0.0011) for set ii1 , (600, 0.0000332) and (300, 0.000132683) for set iii3 , (600, 0.00001741) and (300, 0.0000695824) for set iii2 , and (600, 0.00000693) and (300, 0.0000277) for set iii1 . Figs. 1 and 2 plot P (s; r , L , f ) and sτ P (s; r , L , f ) as a function of s and s/ L D , respectively, where τ = 1.27 and D = 2.75 [13]. The 18 curves in Fig. 1 almost fall on the nine collapsed curves shown in Fig. 2 where each coalesced curve corresponds to a value of L f ν and a fixed r. These simulations suggest the existence of the scaling function G. Assume that the characteristic length of a system with aspect ratio r satisfies
ξ(r , f ) = A r f −ν ,
227
(2)
where A r is a nonuniversal metric factor for a given r. Consider s = sP (s; r , L , f )ds and simulate two cases: (a) f = 0 and a given finite L (i.e., L ξ ), which corresponds to boundary dissipation dominating the system. In this case, L determines system behavior. (b) L → ∞ and a given finite f > 0 (i.e., ξ L), which corresponds to bulk dissipation dominating the system. In this case, ξ determines system behavior. For case (a), consider r = 1, 2, and 3. Simulating P (s; r , L , f = 0) for L = 300, 600, 900, and 1200 shows that s = cr L σ . For case (b), choose a large L to represent the condition L → ∞. Then, select L = 1200 and set r = 1, 2, and 3. Simulating P (s; r , L , f ) for f = 10−1 , 10−2 , 10−3 , and 10−4 shows that s = er f −β . Taken together, cases (a) and (b) imply that
s =
cr L σ
dr ξ σ = er f −β
for L ξ,
(3)
for ξ L .
Previous research shows σ = 1.97 [13]. The knowledge of critical phenomena suggests that σ = σ . Then, β = σ ν . Numerical simulations give that the terms cr , dr , and er are constants for a fixed r while er = e is independent of r. Assuming that σ σ
cr
c1
=
dr
d1
and
dr
d1
=
er /[ A r ]
e 1 /[ A 1 ]σ
=
A1 Ar
(4)
leads to the following relationship
Ar =
c1 cr
1/σ A1.
(5)
Since the value of A r / A 1 is important and the real vale of A 1 is irrelevant in this study, set A 1 = 1 for simplification leading to A r = [c 1 /cr ]1/σ . Thus, this study employs the values of cr obtained from the first moment analysis (i.e., s ) and Eq. (5) to calculate the values of A r . This leads to values of (cr , A r ) = (0.114, 1.000) for r = 1, (0.046, 1.582) for r = 2, and (0.024, 2.199) for r = 3 where the error bars are within 0.005. To unify G with different r into a universal scaling function, this study proposes
P (s; r , L , f ) = s−τ H s/[ L / A r ] D , [ L / A r ]/ f −ν ,
(6)
where H is a two-variable function. To verify Eq. (6), Fig. 3 plots sτ P as a function of s/( L / A r ) D . The 18 curves in Fig. 1 almost fall on the 3 coalesced curves shown in Fig. 3 where sets ir , iir , and iiir for r = 1, 2, and 3 fall on sets I, II, and III, respectively. Note that sets I, II, and III correspond to L /ξ = [ L / A r ]/ f −ν ) = 50, 10, and 1.58, respectively. This suggests the existence of the universal scaling function H . This study uses the critical phenomena framework to investigate a SOC model with different aspect ratios, where a losing probability f functions as a tuning parameter and f = 0 simulates the critical point of the system. Based on the connection between f and ξ , this study proposes the universal scaling function form. The proposed procedure is applicable to calculating the universal scaling functions of other physical quantities of sandpiles, e.g., the toppling area or duration of an avalanche. Finally, fitting exponents τ and D requires the preparations of several (L, f , A r ) which satisfy the condition that [ L / A r ] f ν is a constant according to the expression of Eq. (6). Otherwise, the obtained exponent values will be inaccurate.
228
C.-F. Chen et al. / Computer Physics Communications 182 (2011) 226–228
References
Fig. 3. sτ P (s; r , L , f ) as a function of s/[ L / A r ] D .
Acknowledgement C.-Y. Lin acknowledges the support of the National Science Council of Taiwan, R.O.C. under Grant No. NSC 97-2112-M-194003-MY2.
[1] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381; P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. A 38 (1988) 364. [2] H.J. Jensen, Self-Organized Criticality, Cambridge University Press, New York, 1998. [3] D. Sornette, A. Johansen, I. Dornic, J. Phys. I 5 (1995) 325. [4] A. Vespignani, S. Zapperi, V. Loreto, Phys. Rev. Lett. 77 (1996) 4560. [5] A. Vespignani, S. Zapperi, Phys. Rev. Lett. 78 (1997) 4793; A. Vespignani, S. Zapperi, Phys. Rev. E 57 (1998) 6345. [6] R. Dickman, A. Vespignani, S. Zapperi, Phys. Rev. E 57 (1998) 5095; A. Vespignani, R. Dickman, M.A. Muñoz, S. Zapperi, Phys. Rev. Lett. 81 (1998) 5676. [7] M. Paczuski, S. Boettcher, Phys. Rev. Lett. 77 (1996) 111. [8] L. Pietronero, A. Vespignani, S. Zapperi, Phys. Rev. Lett. 72 (1994) 1690. [9] E.V. Ivashkevich, Phys. Rev. Lett. 76 (1996) 3368. [10] C.-Y. Lin, A.-C. Cheng, T.-M. Liaw, Phys. Rev. E 76 (2007) 041114; A.-C. Cheng, C.-F. Chen, C.-Y. Lin, Int. J. Mod. Phys. C 19 (2008) 1695; C.-Y. Lin, Phys. Rev. E 81 (2010) 021112. [11] S.S. Manna, J. Phys. A 24 (1991) L363. [12] S. Lübeck, Phys. Rev. E 61 (2000) 204; S. Lübeck, Int. J. Mod. Phys. B 18 (2004) 3977. [13] A. Chessa, A. Vespignani, S. Zapperi, Comput. Phys. Commun. 121–122 (1999) 299. [14] C. Tebaldi, M. De Menech, A.L. Stella, Phys. Rev. Lett. 83 (1999) 3952. [15] Y. Shilo, O. Biham, Phys. Rev. E 67 (2003) 066102. [16] O. Malcai, Y. Shilo, O. Biham, Phys. Rev. E 73 (2006) 056125. [17] C.-Y. Lin, C.-F. Chen, C.-N. Chen, C.-S. Yang, I.-Min Jiang, Phys. Rev. E 74 (2006) 031304; C.-F. Chen, C.-Y. Lin, Int. J. Mod. Phys. C 20 (2009) 273. [18] H.E. Stanley, Introduction to Phase Transitions and Critical Phenomena, Oxford Univ. Press, New York, 1971.