Table Of Content

Sandpile avalanche dynamics on scale-free networks 4 0 0 2 D.-S. Lee, K.-I. Goh, B. Kahng and D. Kim n a J SchoolofPhysicsandCenterforTheoreticalPhysics, SeoulNationalUniversity NS50, 7 Seoul151-747, Korea 2 ] h c e Abstract m - t Avalanche dynamics is an indispensable feature of complex systems. Here we study the a t self-organizedcriticaldynamicsofavalanchesonscale-freenetworkswithdegreeexponent s . g throughtheBak-Tang-Wiesenfeld (BTW)sandpilemodel.Thethresholdheightofanode at iissetask1−h with0≤h <1,wherek isthedegreeofnodei.Usingthebranchingprocess m i i approach,weobtaintheavalanchesizeandthedurationdistributionofsandtoppling,which d- followpower-lawswithexponentst andd ,respectively.Theyaregivenast =(g −2h )/(g − n 1−h )andd =(g −1−h )/(g −2)forg <3−h ,3/2and2forg >3−h ,respectively. The o power-lawdistributions aremodifiedbyalogarithmiccorrection atg =3−h . c [ 1 Keywords: avalanche, scale-free network,branching process v PACS:89.70.+c, 89.75.-k,05.10.-a 1 3 5 1 0 4 1 Introduction 0 / t a m Frequently, complex systems in nature as well as in human society suffer massive - d catastrophes triggered from only asmallfraction oftheirconstituents.Unexpected n epidemic spread of diseases and the power outage in the eastern US of the last o c year are the examples of such avalanche phenomena. Such a cascading dynam- : v ics is not always harmful to us. The information cascades making popular hits of i books, movies, and albums are good to writers, actors, and singers, respectively. X Thus it is interesting to understand and predict how those cascades propagate in r a complex system. Recently, the network approach, by which a system is viewed as a network consisting of nodes representing its constituents and links interactions between them,simplifiescomplicateddetailsofcomplexsystems.Such asimplifi- cationunveilsahiddenordersuchasscale-freebehaviorinthedegreedistribution. Here degree is the number of links connected to a certain node. The Internet at theautonomoussystemlevel,theWorld-WideWeb, socialacquaintance networks, biological networks, and other many complex networks exhibit power-law degree PreprintsubmittedtoElsevierScience 2February2008 distributions, p (k)∼ k−g . The networks following such power-law degree distri- d butions are called scale-free (SF) networks [1], where non-negligible fractions of hubs,thenodeswithextraneouslylargedegrees,exist. In this paper, we investigatethe avalanche dynamics on such SF networks through the Bak-Tang-Wiesenfeld (BTW) sandpile model [2], a prototypical model ex- hibitingself-organized criticality(SOC). The studyof sandpiledynamicshas been carriedoutmostlyonregularlatticesintheEuclideanspace.Inthestationarystate, whichcanbereachedwithouttuningaparameter,thesystemexhibitscale-invariant features in the power-law form of the avalanchesize distribution p (s) and the du- a rationorlifetimedistributionℓ(t)as p (s)∼s−t and ℓ(t)∼t−d . (1) a Recently, Bonabeau has studied the sandpile dynamics on the Erdo˝s-Reńyi (ER) randomnetworks[3]andfoundthattheavalanchesizedistributionfollowsapower lawwiththeexponentt ≃1.5,consistentwiththemean-fieldsolution[4].Recently, Lise and Paczuski [5] studied the Olami-Feder-Christensen model [6] on regular ER networks, where degree of each node is uniform but connections are random. Theyfoundtheexponenttobet ≈1.65.However,whendegreeofeachnodeisnot uniform, they found no criticality in theavalanche size distribution.Note that they assumedthatthethresholdofeachnodeisuniform,whereasdegreeisnot.Herewe studytheBTWsandpilemodelon SF networks,wherethethresholdz ofthenode i i is givenas k1−h withk thedegreeofi and 0≤h <1. Wefind thattheexponents i i for theavalanche size and the duration distributiondepend on thedegree exponent g as t =(g −2h )/(g −1−h ) and d =(g −1−h )/(g −2) for g <3−h while, for g >3−h , they show the same behaviors as the conventional mean-field solutions as observedfortheERrandomnetworks. 2 Sandpile model Wepresent thedynamicruleoftheBTWsandpilemodelon agivennetwork. (1) Each nodei is givenaprescribed threshold z (≤k). Thesmallestintegernot i i smallerthan z is denotedas ⌈z ⌉(⌈z ⌉≤k). i i i i (2) At each time step, a grain is added at a randomly chosen node i. The integer- valuedheightofthenodei, h, increases by 1. i (3) If the height at the node i reaches or exceeds z , then it becomes unstableand i the ⌈z ⌉ grains at the node topple to its ⌈z ⌉ randomly chosen adjacent nodes i i amongk ones; i h →h −⌈z ⌉, and h =h +1 forall nodes j which arechosen. i i i j j (4) If this topplingcauses any of the adjacent nodes receiving grains to be unsta- ble, subsequent topplings follow on those nodes in parallel until there is no 2 unstablenodeleft.Thisprocess defines an avalanche. (5) Repeat (2)–(4). Here thethresholdz ofnodei isgivenas i z =k1−h (0≤h <1), (2) i i which is ageneralization ofz =k previouslyinvestigatedin Ref. [7]. Weconcen- i i trate on the distributions of (i) the avalanche area A, i.e., the number of distinct nodes involving in a given avalanche, (ii) the avalanche size S, i.e., the number of topplingeventsinagivenavalanche, and (iii)thedurationT ofagivenavalanche. 3 Branching process approach The mapping of each avalanche to a tree provides a useful way of understanding thestatisticsofavalanchedynamicsanalytically.Foreachavalancheevent,onecan draw a corresponding tree: The node where an avalanche is triggered corresponds totheoriginatorofthetreeandthefollowingnodestodescendants.Inthetreestruc- ture, a descendant born at time t is located away from the originator by distance t along the shortest pathway. The tree stops to grow when no further avalanche pro- ceeds. Then the ensemble of avalanches can be identified with that of trees grown through the branching process. In this mapping, the avalanche duration T is equal to the lifetime of the tree minus one, and the avalanche size S differs from the tree size only by the number of boundary nodes of the tree, which is relatively small whentheoveralltreesizeisverylarge.Ifoneassumesthatbranchingeventsatdif- ferent nodes occur independently and that there is no loop in the tree, the tree size and lifetimedistributioncan be obtained analytically [8,9]. Thosedistributionsare expected to share the same asymptotic behaviors with the avalanche size and duration distribution,respectively, due to the near-equivalence between an avalanche and itscorrespondingtreeintheirscales as mentionedabove. In the branching process describing an avalanche, after initial branching into k de- scendantswithprobabilityq (k),successivebranchingsareassumedtooccurinde- 0 pendentlywithprobabilityq(k).q (k)andq(k)maybedifferentingeneral,butthe 0 statisticsoftheoverallsizeanddurationofanavalancheisdetermineddominantly byq(k).Wecheckedalsonumericallythecasewhereanewgrainisaddedtoanode with the probability proportional to the degree of that node, which gives different q (k)fromthatinthecasewhereanewgrainisaddedrandomly,andfoundthatthe 0 nature of the avalanche dynamics is the same in both cases. Thus, for simplicity, we consider the branching process where every branching occur with probability q(k).FortheBTWmodelintheEuclideanspace,wherethethresholdz ofnodeiis i equaltoitsdegreek,q(k)hasafinitecut-offsuchthatq(k)=0fork>z =const, i i becausethedegreeofeachnodeisuniformandfinite.Consequently,theexponents 3 of the avalanche size and the duration distributions in Eq. (1) come out to be the so-called mean-field values; t = 3/2 and d = 2 [8,9]. These results are known to holdfortheBTWmodelonregularlatticeswithdimensionslargerthan4[4].Note that when dimension is smaller than 4, the branching process approach cannot be applied,so thatthevaluesoftheexponentst and d wouldnotbetrivial. In SF networks, avalanches usually do not form loops, generating tree-structures: According to the numerical simulations of the BTW model for the case of z = k i i on SF networks [7], the statistics of the two quantities A and S are nearly equal when theyare large: Forexample,themaximumarea and size(A , S ) among max max avalanches are (5127, 5128), (12058, 12059) and (19692, 19692) for scale-free networks with g =2.01, 3.0, and ¥ , respectively. The fact that A and S are almost thesameimpliesthattheavalanchestructurecanbetreatedasatree.Fromnowon, weshallnotdistinguishAandS,denotedbys.Thusitisvalidtousethebranching process approach tounderstandtheavalanchedynamicson SF networks. We study the BTW model on SF networks with the degree exponent g and the threshold given as Eq. (2). The branching probability q(k) consists of two factors, that is, q(k) = q (k)q (k), where q (k) is the probability that the threshold z of 1 2 1 i node i is in the range k−1 < z ≤ k and q (k) is the probability that the total i 2 number of grains at the node reaches or exceeds the threshold. If z = f(k) with i i f(x) a monotonic increasing function of x satisfying f(x) ≤ x for all x ≥ 1, the condition of k−1 < z ≤ k implies that q (k) is nothing but the probability that i 1 a node i connected to the one end of a randomly chosen edge has its degree k i in the region(f−1(k−1), f−1(k)],and thus q (k)=(cid:229) ⌊f−1(k)⌋ k′p (k′)/hki, 1 k′=⌊f−1(k−1)⌋+1 d where ⌊x⌋ is the largest integer not larger than x. Notice that (cid:229) ¥ q (k) = 1 and k=1 1 q (k)∼k(1−g +h )/(1−h ) for large k if f(x)≃x1−h (0≤h <1) for large x. q (k) is 1 2 the probability that the node i has height k−1 at the moment of receiving a grain from one of its neighbors. We have checked numerically that a typical height of node is absent, so that all possible k values 0,1,...,k−1 are equally likely [10]. Thus we set q (k)=1/k. As a result, the branching probability q(k) for large k is 2 givenasymptoticallyas 1 g −2h q(k)= q (k)∼k−g ′ g ′ = . (3) k 1 (cid:18) 1−h (cid:19) Whenz =k orh =0,g ′ isreducedtog .Sinceweareinterestedinthecaseofg >2 i i and 0≤h <1, g ′ >0. Usingtheindependenceofthebranchingsfromdifferentparent-nodes,onecande- rivethefollowingself-consistentrelationforthetreesizedistribution p(s) as[8,9] ¥ ¥ ¥ ¥ p(s)=q(0)d + (cid:229) q(k) (cid:229) (cid:229) ··· (cid:229) p(s )p(s )...p(s )d . (4) s,1 1 2 k (cid:229) ki=1si,s−1 k=1 s1=1s2=1 sk=1 4 This relation can bewritten in a morecompact form by introducingthegenerating functions,P(y)=(cid:229) ¥ p(s)ys and Q (w )=(cid:229) ¥ q(k)w k as s=1 k=0 P(y)=y Q (P(y)). (5) Then w =P(y) isobtainedby invertingy=P−1(w )=w /Q (w ). The average size hsi of a finite tree can be obtained easily from the generating functions. finite (cid:229) hsi= sp(s)=P′(1), (6) s=1 whereP′(y)=dP(y)/dy. Usingtherelation,Eq. (5), weobtain Q (P(1)) hsi=P′(1)= , (7) 1−Q ′(P(1)) whereagain Q ′(w )=dQ (w )/dw . The distribution of duration, i.e., the lifetime of the tree can be evaluated simi- larly [8,9]. Let r(t) be the probability that a branching process stops at or prior to timet.ThenfollowingthesimilarstepsleadingtoEq.(4),i.e.,r(t)=(cid:229) ¥ q [r(t− k=0 k 1)]k, onehas r(t)=Q (r(t−1)). (8) Forlarget,r(t)comescloseto1.Onecanobtainw =r(t−1)bysolvingdw /dt ≃ r(t)−r(t−1)=Q(w )−w . Then the lifetime distribution ℓ(t) is obtained through ℓ(t)=r(t)−r(t−1)≃dw /dt. 4 Avalanche sizeandduration distribution Thegrowthofa treedependson theaverage numberofbranches defined as ¥ (cid:229) C= kq(k). (9) k=1 When C > 1 (C < 1), a tree can (cannot) grow infinitely in a probabilistic sense. Thus the case of C = 1 is a critical point for the growth of a tree. One can see thatforanybranchingprocess withq(k)=(1/k)q (k)(k≥1)and(cid:229) ¥ q (k)=1, 1 k=1 1 the average number of branches C is always 1, independent of detailed structural properties of networks. Therefore our assumption q (k)= 1/k corresponds to the 2 conditionfortheself-organized criticality(SOC) ofthesandpilemodel. The inverse function P−1(w ) satisfies P−1(1) = 1. When C = 1, the first-order derivative¶ P−1(w )/w¶ at w =1 is zero and thus P(y) becomes singular at y=1. P(y)isexpandedaroundy=1asP(y)≃1−b(1−y)f withconstantband0<f < 5 1.Thentheasymptoticbehavioroftheavalanchesizedistribution p(s)forlargesis givenby p(s)∼s−f −1,becauseifaseries(cid:229) ¥ a ys withtheradiusofconvergence s=0 s 1 hastheasymptoticbehavior ¥ (cid:229) a ys ∼(1−y)f as y→1, then a ∼s−f −1 as s→¥ . (10) s s s=0 Thefunctionalformofthebranchingprobabilityq(k)determinesthesingularityof P(y). To illustratethis,wefirst considerasimplecasethat 1−a (k=0),  q(k)= a (k=2), (11)    0 (otherwise),     where 0 < a < 1. Then the average number of branches C = (cid:229) kq(k) = 2a and thegeneratingfunctionQ (w )=(cid:229) ¥ q(k)w k =1−a+aw 2.Usingtherelationsof k=0 y=w /Q (w )and w =P(y), itis obtainedthat 1− 1−4a(1−a)y2 P(y)= . (12) p 2ay ThevalueofP(1)=(cid:229) finitep(s)is givenas s=1 1−|1−2a| 1 for 0<a≤ 1 (C≤1), P(1)= = 2 (13) 2a  1−a for 1 <a<1 (C>1), a 2  which means that when 1/2<a<1 (C>1), a tree can grow infinitely with probability 1−P(1) = (2a−1)/a, and the critical point is located at a = 1/2. Near c y=1,P(y)≈1− 2(1−y)fromEq.(12),leadingtof =1/2.Then,theavalanche sizedistribution pp(s) behavesas p(s)∼s−3/2. Ontheotherhand,usingEq.(7), 1 (a<a ), hsi=P′(1)= 2(ac−a) c (14)  1−a (a>a ). 2a(a−ac) c  Even for the case that q(k) has a finite cut-off larger than 2 or decays exponen- tially, the above result holds. This is the conventional mean-field solution for the avalanchesize distribution[4,8,9] and has been shownto holdfor theBTW model on theERrandom networks[3]. When q(k) decays slowly as in Eq. (3), however, its generating function Q (w ) is singularatw =1.Forq(k)inEq.(3),theexpansionofQ (w )aroundw =1isgiven 6 as A (1−w )g ′−1 (2<g <g ), 1 c  Q (w )≃1−(1−w )+ −A2(1−w )2ln(1−w ) (g =g c), (15)  A (1−w )2 (g >g ), 3 c     where A ’s are constants, g ′ is given in Eq. (3), and g = 3−h . The derivation of i c thelogarithmiccorrectionforthecaseofg =g canbefoundin [11].Notethatthe c singularterm (1−w )g ′−1 is thesecondleading term of1−Q (w ) for g <g . Using c the relation P−1(w ) = w /Q (w ) in Eq. (5), the behavior of P(y) around y = 1 is obtained for each region of g from Eq. (15), and in turn, using Eq. (10), p(s) for s→¥ . Wefind that s−(g −2h )/(g −1−h ) (2<g <g ), c   p(s)∼s−3/2(lns)−1/2 (g =g ), (16)  c  s−3/2 (g >g c).     Thus,theexponentt isgivenast =(g −2h )/(g −1−h )for2<g <g andt =3/2 c forg ≥g . c Alsoobtainedisr(t)fromEq.(15)byusingEq.(8).Thedurationdistributionℓ(t), whichis thederivativeofr(t),isfound tobe t−(g −1−h )/(g −2) (2<g <g ), c   ℓ(t)∼t−2(lnt)−1 (g =g ), (17)  c  t−2 (g >g c).     Thatis,theexponentd isgivenasd =(g −1−h )/(g −2) for2<g <g andd =2 c forg ≥g . c 5 Conclusion WehavestudiedtheBTWsandpilemodelonSFnetworkswiththedegreeexponent g to understand the avalanche dynamics in complex systems. The main results are the avalanche size and duration distribution. The exponents t and d increase with increasing g , implying that the hubs play a role of reservoir, that is, sustain large amountofgrainstomaketheSFnetworkresilientunderavalanchedynamics.This isreminiscentofthestructuralresilienceoftheSF networkunderrandomremoval 7 of nodes for g ≤ 3 [12,13,14]. We also checked the case where the threshold z i containsnoiseinthewaythatz =z k withz beingdistributeduniformlyin[0,1]. i i i i Wefindthatsuchavariationdoesnotchangethenatureoftheavalanchedynamics. However, when the threshold is given in terms of a quantity other than degree, e.g., load,thecorrespondingavalanchedynamicshas no reason tofollowthesame statisticsas studiedinthispaper, which remainsfurtherworks. This work is supported by the KOSEF Grant No. R14-2002-059-01000-0 in the ABRL program. References [1] A.-L.Baraba´siandR.Albert,Science286,509(1999). [2] P.Bak,C.Tang,andK.Wiesenfeld, Phys.Rev.Lett.59,381(1987);Phys.Rev.A38, 364(1988). [3] E.Bonabeau, J.Phys.Soc.Japan64,327(1995). [4] P.Alstrøm,Phys.Rev.A38,4905(1988). [5] S.LiseandM.Paczuski,Phys.Rev.Lett.88,228301(2002). [6] Z. Olami, H.J.S. Feder, and K. Christensen, Phys. Rev. Lett. 68, 1244 (1992); K. Christensen andZ.Olami,Phys.Rev.A46,1829(1992). [7] K.-I.Goh,D.-S.Lee,B.Kahng,andD.Kim,Phys.Rev.Lett.91,148701(2003). [8] T.E.Harris,TheTheoryofBranchingProcesses(Springer-Verlag, Berlin,1963). [9] R.Otter,Ann.Math.Statist.20,206(1949). [10] D.-S.Lee,K.-I.Goh,B.Kahng,andD.Kim,toappearinJ.KoreanPhys.Soc.(2004). [11] J.E.Robinson, Phys.Rev.83,678(1951). [12] R.Albert,H.Jeong, andA.-L.Baraba´si, Nature(London)406,378(2000). [13] R.Cohen,K.Erez,D.ben-Avraham, andS.Havlin,Phys.Rev.Lett.85,4626(2000). [14] D.S. Callaway, M.E.J. Newman, S.H. Strogatz, and D.J. Watts, Phys. Rev. Lett. 85, 5468(2000). 8