Critical behavior and conservation in directed sandpiles

Critical behavior and conservation in directed sandpiles Romualdo Pastor-Satorras and Alessandro Vespignani The Abdus Salam International Centre for Theoretical Physics, P.O. Box 586, 34100 Trieste, Italy 共Received 20 June 2000兲 We perform large-scale simulations of directed sandpile models with both deterministic and stochastic toppling rules. Our results show the existence of two distinct universality classes. We also provide numerical simulations of directed models in the presence of bulk dissipation. The numerical results indicate that the way in which dissipation is implemented is irrelevant for the determination of the critical behavior. The analysis of the self-affine properties of avalanches shows the existence of a subset of superuniversal exponents, whose value is independent of the universality class. This feature is accounted for by means of a phenomenological description of the energy balance condition in these models. PACS number共s兲: 05.65.⫹b, 05.70.Ln I. INTRODUCTION Sandpile cellular automata are the most famous example of self-organized critical 共SOC兲behavior 关1–3兴. Under an external drive consisting of a slow addition of sand 共energy兲 grains and the action of dissipation through the loss of energy on the lattice boundaries, these models reach a stationary steady state. In the limit of infinitesimal driving and dissipation 共this last achieved in the thermodynamic limit兲, the stationary state of sandpile models exhibits diverging response functions associated to a characteristic avalanche dynamics. This is the hallmark of a critical behavior that has attracted an enormous amount of interest as a plausible explanation of the avalanche-like critical behavior empirically observed in many natural systems 关3兴. Sandpile models have been at the center of an intense research activity made of both analytical studies and numerical simulations. Despite the simple definition of these automata, it turns out that their full analytical understanding is a very problematic task 关4兴. As a further complication, also the numerical inspection of these models results to be particularly difficult. For example, the precise identification of universality classes has resisted for many years even the most careful numerical analysis, and only recent results have partially settled this problem 关5–8兴. On the other hand, these refined analyses have pointed out that several sandpile models do not follow the simple finite size scaling 共FSS兲form usually adopted in the description of critical behavior 关9兴. For instance, the more sophisticated multiscaling approach 关10–12兴seems to be required for a full description of the scaling properties of the original Bak, Tang, and Wiesenfeld 共BTW兲model 关1,2兴. Many sandpile features have been underlined as the possible origin of these scaling anomalies. The deterministic dynamical rules of the BTW model induce nonergodic effects 关8兴, that are certainly missing in stochastic models, such as the Manna model 关13,4兴, which shows a perfect FSS behavior, even for moderate system sizes. A further complication of sandpile automata stems from the peculiar role of the boundary dissipation, that makes the lattice size scaling entangled with the system dynamics. In such cases, the thermodynamic limit is essential for the dissipative dynamics of large avalanches. A clear understanding of the interplay between dissipation and size scaling has not yet been achieved and it has been recently the subject of several studies 关11,14兴. In this paper we address some of the aforementioned problems in the case of directed sandpile models 关10,15–18兴. In this case Dhar and Ramaswamy obtained an exact solution for the Abelian deterministic directed sandpile 共DDS兲关15兴, that can be used as a touchstone to check the numerical simulation analysis. Directed sandpiles thus become an interesting test field to study how the critical behavior is affected by the introduction of stochastic elements and dissipation. We perform large scale numerical simulations of two directed sandpile automata: the deterministic directed sandpile model 关15兴and the stochastic directed sandpile model 关18兴. We study both models in the case of boundary and bulk dissipation 关19–22兴. We find, in agreement with the results in Ref. 关18兴, that the models define two different universality classes. In addition we show that the universality class of the models does not depend on the way in which dissipation is implemented. Finally we analyze the properties of anisotropic models in which the dynamics is not fully directed 关10,23兴. In this case we observe that on large scales the critical behavior is the same of that of fully directed models. Results for the stochastic models are compared with a recent theoretical approach by Paczuski and Bassler 关24兴, that provides values for the critical exponents in perfect agreement with numerical simulations. These results are also recovered in Ref. 关25兴. The numerical analysis also points out that some critical exponent values, such as the correlation length exponents or the affinity exponent 共to be defined later on兲, are independent of the particular universality classes and common to all models considered. In order to explain this numerical evidence, we provide a phenomenological characterization of directed sandpiles based on the basic symmetries introduced by the conserved dynamics of these automata. Following balance of energy arguments inspired in Refs. 关26–28兴, we derive a series of results and predictions on the value of critical exponents which are a straightforward consequence of conservation. These general results can be considered as superuniversal, because they characterize the critical behavior of all directed sandpiles with local dynamical rules, independently on the specific universality class. The results presented here provide a general picture of directed models and the role of PHYSICAL REVIEW E NOVEMBER 2000VOLUME 62, NUMBER 5 PRE 62 1063-651X/2000/62共5兲/6195共11兲/$15.00 6195 ©2000 The American Physical Society boundary and bulk dissipation in the process of selforganization. The paper is arranged as follows. In Sec. II we introduce and define the various directed models considered. Sections III and IV present and discuss, from the standpoint of universality, the numerical results for directed models with boundary and bulk dissipation. In Sec. V we introduce anisotropic models, and present the numerical results obtained, in comparison with those of directed models. Section VI is devoted to an analytical approach based on the conservation of energy. Finally, in Sec. VII we draw our conclusions and perspectives. II. DIRECTED MODELS Sandpile models are usually defined on a d-dimensional hypercubic lattice of size L. To each node of the lattice is assigned an integer variable zi, called ‘‘energy.’’ Energy is added to the system uniformly at randomly chosen sites (zi →zi⫹1). When a site becomes active, that is, when its energy becomes larger than or equal to a certain threshold zc, it topples. A toppling site loses an energy zc, that is distributed among its neighbors according to a certain set of rules. The neighbors that receive energy can become active and topple on their turn, thus generating an avalanche. The slow driving condition is effectively imposed by stopping the random energy addition during the avalanche spreading. This means that the driving time scale is infinitely large with respect to the toppling characteristic time scale. The models we consider in this section are directed, in the sense that the energy is always transported along a preferred fixed direction. We denote this preferred direction by the coordinate x 储 , whose positive direction is usually defined as ‘‘downwards.’’ The transverse direction 共subspace of dimension d⫺1 perpendicular to x 储 ) will be denoted by x ជ ⬜. The toppling rules of the models define two main classes. 共i兲Deterministic directed sandpile 共DDS兲:Inddimensions, the threshold is set to zc⫽2d⫺1. When a site in a given hyperplane x 储 topples, it sends deterministically one grain of energy to each one of its nearest and next-nearest neighbors on the hyperplane x 储 ⫹1关see Fig. 1共a兲兴. Our definition is somewhat different from the original model of Dhar and Ramaswamy 关15兴, in both the driving and the orientation of the lattice. Both models, however, are expected to share the same universality class, being deterministic and directed. Numerical simulations confirm indeed this point 关18兴. 共ii兲Stochastic directed sandpile 共SDS兲: In this case, the threshold is zc⫽2, independently of the dimensionality of the lattice. When a site in the hyperplane x 储 topples, it sends two grains of energy to two sites, randomly chosen among its 2d⫺1 nearest and next-nearest neighbors on the hyperplane x 储 ⫹1. The toppling rules of this model can be defined exclusive if the two energy grains are always distributed on different sites, Fig. 1共b兲. On the other hand, the model can be defined nonexclusive if the dynamics allows the transfer of two energy grains onto the same site, Fig. 1共c兲. We therefore report simulations on the exclusive stochastic directed sandpile 共ESDS兲and on the nonexclusive stochastic directed sandpile 共NESDS兲. In spite of the stochastic nature of these models, we must bear in mind that they are nevertheless Abelian 关4兴. The discussion therefore focuses on the difference between stochastic and deterministic models. Once the toppling rules have been determined, the models are finally defined by specifying the dissipation mechanism. For systems with boundary dissipation, we impose periodic boundary conditions in the transverse directions x ជ ⬜and open at the hyperplane x 储 ⫽L. In this way, the models are locally conserved; energy can only leave the system at the bottom of the lattice. In models with bulk dissipation, we impose periodic boundary conditions in both the x 储 and x ជ ⬜directions. Dissipation is implemented by allowing a toppling site to lose an energy zcwithout transferring it with probability p 关26,21兴. This means that, on average, an energy ⑀ ⫽zcpis dissipated in each toppling. In the limit ⑀ →0, the system shows critical behavior 关26兴. In the stationary state we can define the probability that the addition of a single energy grain is followed by an avalanche of toppling events. Avalanches are then characterized by the total number of topplings sand the time duration t.In the limit of infinitesimal driving 共slow driving condition兲the system shows scaling behavior and the probability distributions of these quantities follow the finite-size scaling 共FSS兲 forms P共s兲⫽s⫺ ␶ sG共s/sc兲,共1兲 P共t兲⫽t⫺ ␶ tF共t/tc兲,共2兲 where scand tcare the characteristic size and time, respectively. The exponents ␶ sand ␶ tcharacterize the critical behavior and define the universality classes to which the models belong. In the critical region the characteristic time and size are determined only by the system size Lor the dissipation ⑀ , in the case of boundary and bulk dissipation, respectively. In directed models, the affinity exponent ␨ is of particular importance; it relates the avalanche characteristic lengths in the perpendicular direction, ␰ ⬜, and in the parallel direction, ␰ 储 , through the relation ␰ ⬜⬃ ␰ 储 ␨ . This exponent characterizes the degree of anisotropy due to the preferential direction present in the transport of the energy. In other words, it expresses the self-affine properties in the scaling of avalanches. A general result concerns the average avalanche size 具 s 典 , that also scales linearly with L关10,15,23兴; this result can be exactly obtained by inspecting the conservation symmetry of the model as we shall see in Sec. VI. For the DDS, the exact analytical solution in d⫽2 yields the exponents ␶ s⫽4/3 and ␶ t⫽D⫽3/2 关15兴. The upper critical dimension is found to be dc⫽3, and it is also possible to FIG. 1. Toppling rules in d⫽2 for directed sandpiles. Filled circles represent active 共toppling兲sites; empty circles are stable sites. In the deterministic model 共a兲an active site sends one grain to each of its three neighbors on the next downwards row. In the stochastic models, exclusive 共b兲and nonexclusive 共c兲, one grain is sent to two randomly chosen downwards neighbors. 6196 PRE 62 ROMUALDO PASTOR-SATORRAS AND ALESSANDRO VESPIGNANI find exactly the logarithmic corrections to scaling 关15,29兴. The introduction of stochastic ingredients in the toppling dynamics of directed sandpiles has been studied recently in a model that randomly stores energy on each toppling 关16兴. This model is strictly related to directed percolation and defines a universality class ‘‘per se.’’ It would be very interesting to investigate the possible connection between this model and the stochastic directed one presented in this work. III. NUMERICAL SIMULATIONS WITH BOUNDARY DISSIPATION In this section we report results from computer simulations of deterministic and stochastic directed sandpiles, performed with boundary dissipation. The system sizes considered range from L⫽100 to L⫽6400. The statistical distribution functions have been computed averaging over 107nonzero avalanches. In the case of boundary dissipation, the lattice size Lis the only characteristic length present in the system. Approaching the thermodynamic limit (L→⬁), the avalanche characteristic size and time in Eqs. 共1兲and 共2兲diverge as sc⬃LDand tc⬃Lz, respectively. The exponent Ddefines the fractal dimension of the avalanche cluster and zis the usual dynamic critical exponent. The directed nature of the model introduces a drastic simplification, since it imposes z⫽1. In order to compute the different exponents characterizing the dynamics of the avalanches, we have performed the moment analysis of the distributions, in analogy to the method developed by De Menech et al. 关11,12兴. We define the qth moment of the avalanche size distribution on a lattice of size L as 具 sq 典 L⫽ 兰 dssqP(s). If the FSS hypothesis 共1兲is valid in the asymptotic limit of large s, then the qth moment has the following dependence on system size: 具 sq 典 L⫽LD(q⫹1⫺ ␶ s) 冕 dy y(q⫺ ␶ s)G共y兲⬃L ␴ s(q).共3兲 The exponent ␴ s(q)⫽D(q⫹1⫺ ␶ s) is computed as the slope of the log-log plot of 具 sq 典 Las a function of L. For large enough values of q关i.e., away from the region where the integral in Eq. 共3兲is dominated by its lower cutoff兴, one can compute the fractal dimension Das the slope of ␴ s(q)asa function of q:D⫽ ⳵␴ s(q)/ ⳵ q. On the other hand, since the first moment must scale linearly with L, we have ␴ s(1)⫽1. Once Dis known we can estimate ␶ susing the relation ␴ s(1)⫽D(2⫺ ␶ s)⫽1. Along the same lines we can obtain the moments of the avalanche time distribution. In this case, 具 tq 典 L⬃L ␴ t(q), with ⳵␴ t(q)/ ⳵ q⫽z. Analogous considerations for small qapply also for the time moment analysis. Here, an estimate of the asymptotic convergence of the numerical results is provided by the constraint z⫽1, that must hold for large enough sizes. Then, the ␶ texponent can be found using the scaling relation (2⫺ ␶ t)⫽ ␴ t(1). Once the exponents have been estimated numerically, we can check the accuracy of the moment analysis’ predictions using the FSS hypothesis. If the FSS hypothesis of Eqs. 共1兲,共2兲is correct, then the plots of the distributions, under the rescaling s→s/LDand P(s)→P(s)LD ␶ sand correspondingly t→t/Lzand P(t)→P(t)Lz ␶ t, should collapse onto the same universal function, for different values of L. In Table I we report the exponents found for the DDS, ESDS, and NESDS models in d⫽2. Figure 2 shows the moments ␴ s(q) and ␴ t(q). Figures 3 and 4 plot the FSS data collapse for sizes and times, respectively. The exponents obtained for the DDS are in perfect agreement with the expected analytical results. This fact supports the idea that the system sizes used in the present work allow to recover the correct asymptotic behavior. Results for the ESDS and NESDS are identical within the error bars, pointing out that these two models are in the same universality class. On the other hand, the obtained exponents prove beyond any doubts that deterministic and stochastic directed sandpile models do not belong to the same universality class. We have also directly computed the characteristic lengths in the parallel and transversal directions ␰ 储 and ␰ ⬜as a funcFIG. 2. Plot of 共a兲 ␴ s(q) and 共b兲 ␴ t(q) for the d⫽2 models DDS, ESDS, and NESDS with boundary dissipation. TABLE I. Critical exponents for directed sandpiles with boundary dissipation in d⫽2. DR: Dhar and Ramaswamy’s exact solution; DDS, deterministic directed model; ESDS and NESDS, stochastic directed models. Figures in parentheses denote statistical uncertainties. Model ␶ sD ␶ tz ␨ DR 4/3 3/2 3/2 1 1/2 DDS 1.34(1) 1.51(1) 1.51(1) 1.00(1) 0.50(1) ESDS 1.43(1) 1.74(1) 1.71(3) 0.99(1) 0.51(1) NESDS 1.43(1) 1.75(1) 1.74(4) 0.99(1) 0.51(1) PRE 62 6197CRITICAL BEHAVIOR AND CONSERVATION IN . . . tion of the system size. The anisotropy of the system is reflected in the different definitions of both characteristic lengths. In this sense, we define them with the same spirit as in directed percolation 关30兴. Consider a given avalanche, labeled ␣ , that has started at the site (x 储 (0) ,x ជ ⬜ (0)), and has affected the set of different sites 兵 (x 储 (i),x ជ ⬜ (i)) 其 , for i⫽0•••a⫺1共i.e., it has covered an area a). Let us define the quantities R 储 共 ␣ 兲⫽1 a兺 i⫽1 a⫺1 兩 x 储 (0)⫺x 储 (i) 兩 共4兲 and R⬜ 2共 ␣ 兲⫽1 a兺 i⫽1 a⫺1 共x ជ ⬜ (0)⫺x ជ ⬜ (i)兲2.共5兲 Furthermore, let us define R 储 (a) and R⬜ 2(a) as the averages of the previous quantities, over all avalanches of the same fixed area a. Let P(a) be the probability of observing an avalanche of area a. We define the correlation lengths by ␰ 储 ⫽兺aR 储 共a兲aP共a兲 兺aaP共a兲, ␰ ⬜ 2⫽ 兺aR⬜ 2共a兲aP共a兲 兺aaP共a兲.共6兲 The different definitions 共4兲and 共5兲are obviously due to the different nature of the avalanche spreading in the directions x 储 and x⬜. In the former case, the spreading is isotropic, and thus the second moment of the relative distance distribution is needed to define a meaningful correlation length. In the latter case, on the other hand, the spreading is always in the direction of growing x 储 , and therefore the first moment is sufficient. The system being critical, both correlation lengths should scale with the system size, defining the exponents ␯ 储 and ␯ ⬜ by ␰ 储 ⬃L ␯ 储 , ␰ ⬜⬃L ␯ ⬜.共7兲 The affinity exponent, defined by ␰ ⬜⬃ ␰ 储 ␨ 共8兲 is thus given by ␨ ⫽ ␯ ⬜/ ␯ 储 . We have calculated the correlations lengths in the models DDS, ESDS, and NESDS, given by the definition 共6兲. The results, plotted in Fig. 5, give the following dependence of the correlation lengths with system size for all models: ␰ 储 ⬃L, ␰ ⬜⬃L1/2.共9兲 These relations define the exponents ␯ 储 ⫽1 and ␯ ⬜⫽1/2, and an affinity exponent ␨ ⫽1/2. It is interesting to note that this exponent is independent of the universality class of the model, defining a sort of superuniversal property of directed models. As pointed out in Ref. 关18兴, the stochastic dynamics of SDS models introduces multiple toppling events on the same site, which are by definition absent in the deterministic case. This gives rise to a very different avalanche structure, evenFIG. 3. Data collapse analysis of the integrated avalanche size distribution for the d⫽2 stochastic models with boundary dissipation 共a兲ESDS and 共b兲NESDS. System sizes are L⫽400, 800, 1600, 3200, and 6400. FIG. 4. Data collapse analysis of the integrated avalanche time distribution for the d⫽2 stochastic models with boundary dissipation 共a兲ESDS and 共b兲NESDS. System sizes are L⫽400, 800, 1600, 3200, and 6400. 6198 PRE 62 ROMUALDO PASTOR-SATORRAS AND ALESSANDRO VESPIGNANI tually reflected in the different asymptotic critical behavior. It is worth remarking that the universality class of SDS appears robust to modifications of the stochastic microscopic dynamics as pointed out in Ref. 关31兴, where it is shown that modifications of SDS models with stochastic toppling threshold still belong to the same universality class. Recently, Paczuski and Bassler 关24兴, have proposed a theoretical approach that allows the calculation of critical exponents in directed models with multiple topplings. The analysis goes through the mapping of the avalanche evolution into the dynamics of an interface moving in a random medium, as also proposed in Refs. 关32,33兴. This theoretical result gives the exponents ␶ s⫽10/7 and ␶ t⫽D⫽7/4, in perfect agreement with the values obtained by numerical simulations, Table I. The same exponent values are also found in the approach of Ref. 关25兴. IV. NUMERICAL SIMULATIONS WITH BULK DISSIPATION In this section we report results from computer simulations of deterministic and stochastic sandpiles, performed with bulk dissipation. In this case, dissipation is implemented as described in Sec. II. That is, in a system with periodic boundary conditions, each toppling site has a probability ⑀ /zcof losing an energy zc, and a probability 1 ⫺ ⑀ /zcof transferring it to its neighbors. The dissipation rates range from ⑀ ⫽0.0016 to 0.0512, and the 共fixed兲system size considered is L⫽6400. Statistical distribution functions have been computed averaging over 107nonzero avalanches. In the presence of bulk dissipation the characteristic sizes are determined by the dissipation rate ⑀ , which defines the only characteristic length in the system. Approaching the limit ⑀ →0, the avalanche characteristic size and time diverge as sc⬃ ⑀ ⫺⌬sand tc⬃ ⑀ ⫺⌬t, respectively. It is also very easy to relate the mean avalanche size to the dissipation rate ⑀ . On average, each added grain must be dissipated in the evolution of the avalanche, resulting in ⑀ 具 s 典 ⫽1. This readily yields 具 s 典 ⫽ ⑀ ⫺1. In this case it is extremely important that the characteristic length of the avalanche ␰ 储 is always smaller than the size of the lattice used. This allows us to study only finite size effects introduced by the dissipation probability, without spurious effects due to the finite lattice size. The moment analysis can be straightforwardly generalized to systems with bulk dissipation. In this case the role of the system size Las scaling parameter is played by the dissipation ⑀ . If the FSS hypothesis holds, the qth moment for, say the size distribution, has an explicit dependence on the dissipation rate that reads 具 sq 典 ⑀ ⬃ ⑀ ⫺⌬s(q⫹1⫺ ␶ s)⫽ ⑀ ⫺ ␳ s(q).共10兲 The new moment ␳ s(q)⫽⌬s(q⫹1⫺ ␶ s) can be estimated by linear regression in a log-log plot of 具 sq 典 ⑀ as a function of ⑀ ⫺1. Once this moment is computed, the exponent ⌬sis given by ⌬s⫽ ⳵␳ s(q)/ ⳵ q. The relation 具 s 典 ⫽ ⑀ ⫺1imposes ␳ s(1)⫽1, and from here, once known ⌬s, we compute ␶ s using the relation ␳ s(1)⫽⌬s(2⫺ ␶ s). Analogous considerations allow us to compute the exponents of the time distribution ⌬tand ␶ t. Finally, to check the exponents with the data collapse technique, one must plot the rescaled functions P(s) ⑀ ⫺⌬s ␶ sas a function of s/ ⑀ ⫺⌬sand P(t) ⑀ ⫺⌬t ␶ tas a function of t/ ⑀ ⫺⌬t, respectively. FIG. 5. Correlation lengths ␰ 储 and ␰ ⬜as a function of Lfor the models with boundary dissipation DDS (䊊), ESDS (䉭), and NESDS (〫). The dashed lines are guides to the eye with slope 1.00 and 0.50. FIG. 6. Plot of 共a兲 ␳ s(q) and 共b兲 ␳ t(q) for the d⫽2 models DDS, ESDS, and NESDS with bulk dissipation. TABLE II. Critical exponents for directed sandpiles with bulk dissipation in d⫽2. DR: Dhar and Ramaswamy’s exact solution; DDS, deterministic directed model; ESDS and NESDS, stochastic directed models. Figures in parentheses denote statistical uncertainties. Model ␶ s⌬s ␶ t⌬t ␨ DR 4/3 3/2 3/2 1 1/2 DDS 1.32(1) 1.50(1) 1.52(1) 1.00(1) 0.51(1) ESDS 1.42(1) 1.72(2) 1.70(4) 0.98(2) 0.51(1) NESDS 1.43(1) 1.75(2) 1.70(5) 0.99(2) 0.50(1) PRE 62 6199CRITICAL BEHAVIOR AND CONSERVATION IN . . . In Table II we report the exponents computed in d⫽2 for the directed models DDS, ESDS, and NESDS with bulk dissipation. The corresponding moments ␳ s(q) and ␳ t(q) are shown in Figs. 6, while Figs. 7 and 8 plot the data collapse for sizes and times, respectively. To conclude our analysis of directed sandpiles with bulk dissipation, we have proceeded to compute the correlation length of the models. In this case, the scaling of the correlation lengths with vanishing dissipation define the scaling exponents ␰ 储 ⬃ ⑀ ⫺ ␯ 储 ⬘, ␰ ⬜⬃ ⑀ ⫺ ␯ ⬜ ⬘,共11兲 and an affinity exponent ␨ ⫽ ␯ ⬜ ⬘/ ␯ 储 ⬘. Using an analogous definition as in the case of boundary dissipation, we compute the exponents ␯ 储 ⬘⫽1, ␯ ⬜ ⬘⫽1/2, and ␨ ⫽1/2, as shown in Fig. 9. That is, the correlation length exponents are identical for both boundary and bulk dissipation. These results again imply an affinity exponent ␨ ⫽1/2 in all the models studied so far.These results confirm that the critical behavior of models with boundary or bulk dissipation is identical. In fact, all critical exponents ␶ s, ␶ t,z, and ␨ are equal in both cases 关34兴. This further confirms the complete equivalence of both points of view with respect to sandpiles and shows that, at least in the directed case, the open boundary conditions usually implemented in simulations do not affect the scaling behavior in a peculiar way. Of course, the open boundary conditions breaks the translational invariance of the system, but in the thermodynamic limit this effect is negligible for the asymptotic critical behavior. Finally, these results validate theoretical approaches in which it is assumed a homogeneous dissipation that is much easier to treat analytically. As a last observation it is worth remarking that also in this case, a series of exponents such as ␨ and ␯ ⬜ ⬘assume values independently of the universality class of the model under study. This sort of superuniversality can be explained in terms of energy conservation as we shall see in Sec. VI. V. NUMERICAL SIMULATIONS OF ANISOTROPIC MODELS An important question to study in directed sandpile models is the effect on the scaling properties of any amount of FIG. 7. Data collapse analysis of the integrated avalanche size distribution for the d⫽2 stochastic models with bulk dissipation 共a兲 ESDS and 共b兲NESDS. Dissipations are ⑀ ⫽0.0256, 0.0128, 0.0064, 0.0032, and 0.0016. FIG. 8. Data collapse analysis of the integrated avalanche time distribution for the d⫽2 stochastic models with bulk dissipation 共a兲 ESDS and 共b兲NESDS. Dissipations are ⑀ ⫽0.0256, 0.0128, 0.0064, 0.0032, and 0.0016. FIG. 9. Correlation lengths ␰ 储 and ␰ ⬜as a function of ⑀ for the models with bulk dissipation DDS (䊊), ESDS (䉭), and NESDS (〫). The dashed lines are guides to the eye with slope 1.00 and 0.50. 6200 PRE 62 ROMUALDO PASTOR-SATORRAS AND ALESSANDRO VESPIGNANI diffusion along the preferred direction of transport x 储 . One would expect that the broken symmetry introduced by the preferential direction should prevail on large scales, so that the dynamical scaling in directed and simply anisotropic sandpiles become indistinguishable in the thermodynamic limit. This fact hints towards the possibility of a unique universality class for both directed and anisotropic sandpiles. This universality class is determined uniquely by the lack of symmetry along the x 储 direction, and the presence or absence of stochastic elements in the definition of the models. In order to test this conjecture, we have performed numerical simulations of an anisotropic stochastic sandpile model, defined according to the following rules: on a hypercubic lattice of size L, we consider a model with threshold zc⫽2. When a site topples, it sends two grains of energy to two sites, randomly selected among the 2d⫺1 nearest and next-nearest neighbors on the hyperplane x 储 ⫹1, and the nearest neighbor on the hyperplane x 储 ⫺1, see Fig. 10. The rules in this model are defined non-exclusive, in such a way that the same site can receive the two sand grain expelled by an active site. The model is clearly anisotropic, because the probability to transfer energy in the downwards direction is three times larger than in the upwards direction. It would thus correspond to a nonexclusive stochastic anisotropic sandpile 共NESAS兲. We consider only the case of boundary dissipation, performing simulations for sizes ranging from L⫽100 up to 6400, and averaging over 107nonzero avalanches. In Fig. 11 we plot the correlation length ␰ 储 and ␰ ⬜, measured according to the rules given in Eqs. 共6兲. We confirm the expectation that anisotropic models have the same scaling properties, as regards the scaling of the correlation lengths, as directed models with the same deterministic or stochastic ingredients. We have also measured the exponents ␶ s, ␶ t,D, and zfor this model, using the moment analysis technique. The values found are ␶ s⫽1.43(1), D⫽1.75(1), ␶ t⫽1.72(2), z⫽0.98(2). These results, compared with Tables I and II, show that this anisotropic models belongs to the same universality class of the ESDS and NESDS directed models, confirming the irrelevance of the diffusion along the preferred direction x 储 . VI. THE ROLE OF CONSERVATION IN SANDPILE MODELS We have seen in the preceding sections that a subset of critical exponents characterizing the critical behavior of directed and anisotropic models have an interesting superuniversal property; i.e., they are independent of the universality class of the models. In order to understand this feature we perform a theoretical analysis based on the conservation of energy, that is the basic symmetry in standard sandpile automata. We shall see in the following that the superuniversal character of some critical exponents is dictated by simple energy conservation considerations. The use of this approach also allows us to establish a relation between boundary and bulk dissipation models by introducing an effective dissipation that depends on the system size. The avalanche dynamics in sandpile models is implicitly due to the imposed infinite time scale separation between driving and dissipation 关26,27,35兴. In order to devise a theory that can take into account the symmetry introduced by the energy conservation, one must first regularize the rules of the models in such a way that a single time scale is ruling the dynamics. One way to do so is to introduce a nonzero driving rate, defined as the probability per unit time hof a site to receive a grain of energy 关26,35兴. This driving rate plays the role of an external field and leads to the SOC behavior in the limit h→0⫹. On the other hand, given that the toppling rules are conserved, energy can leave the system only at the boundaries. Boundary dissipation is a natural choice in computer simulations. However, it introduces undesirable complications due to its singular character in a local theory. It is therefore convenient to use an homogeneous effective dissipation ⑀ , defined as the average energy lost in each toppling event. As observed in previous sections, one can define models with periodic boundary conditions and built-in bulk dissipation. When constructing the local theory for models with open boundary conditions, the bulk dissipation ⑀ amounts to an effective parameter that is to be related to the system size L.With all these ingredients, we are ready to formulate conservation of energy as a continuous equation. In sandpiles, we define the order parameter ␳ aas the density of active sites 共i.e., whose height z⭓zc). The only dynamics in the model is obviously due to the field ␳ a(x ជ ,t), which is coupled to the local energy density E(x ជ ,t)共i.e., the local density of FIG. 10. Toppling rules in d⫽2 for an anisotropic sandpile. Filled circles represent active 共toppling兲sites; empty circles are stable sites. An active site sends one grain to two randomly chosen sites selected among the three downwards neighbors and the upward nearest neighbor. FIG. 11. Correlation lengths ␰ 储 and ␰ ⬜as a function of Lfor the model with boundary dissipation NESAS. The dashed lines are guides to the eye with slope 1.00 and 0.50. PRE 62 6201CRITICAL BEHAVIOR AND CONSERVATION IN . . . sand grains兲, which enhances or suppresses the generation of new active sites. A Langevin description for sandpile automata is possible by considering the dynamics of the local order-parameter field ␳ a(x ជ ,t) in a coarse-grained picture, bearing in mind that the energy density E(x ជ ,t)isaconserved field. In Refs. 关27,28兴, in analogy with absorbing-state phase transitions 关36,37兴, a pair of coupled dynamical equations for the fields ␳ a(x ជ ,t) and E(x ជ ,t) were proposed. In the following we elucidate the consequences of energy conservation and we focus only on the latter equation. The interested reader can find the full set of equations in Ref. 关28兴.Inthe next subsections we shall consider separately directed and anisotropic models. A. Directed sandpiles We seek a continuous equation for the coarse-grained local density of energy E(x ជ ,t). In the limit of zero driving and dissipation, energy is conserved. Therefore, the evolution equation fulfilled by the local field Eis ⳵ E共x ជ ,t兲 ⳵ t⫽⫺ⵜ ជ •J ជ E⫺ ⑀ ␳ a共x ជ ,t兲⫹h共x ជ ,t兲⫹ ␩ E共x ជ ,t兲.共12兲 The first term simply represents the diffusion of energy; the second term accounts for the dissipation that is associated with every toppling event; the third term represents the external driving. Finally, the last term is a source of stochastic noise, that accounts for the randomness in the flow of energy. The noise term can be generated by the toppling rules in a stochastic model, or by the initial conditions plus the random driving in a deterministic model. We will require the noise to have zero average 具 ␩ E共x ជ ,t兲 典 ⫽0. 共13兲 The noise correlator 具 ␩ E(x ជ ,t) ␩ E(x ជ ⬘,t⬘ 典 is of fundamental importance for the determination of universality classes and the critical behavior of the order parameter. However, for our present purposes we do not need precise knowledge of its analytical form 共for a detailed discussion, see Refs. 关27,28兴兲. The current can be constructed by appealing to the symmetries of the model. The transport of energy is due to topplings. These are isotropic along the transversal direction x ជ ⬜, therefore the current along this direction will be proportional to the gradient of the density of active sites. In the preferred direction, on the other hand, all the energy is transferred downwards; therefore, the current in this direction must be proportional to the density of active sites. The final form of the current is then J ជ E共x ជ ,t兲⫽⫺D⬜ⵜ ជ ⬜ ␳ a共x ជ ,t兲⫹2␭ ␳ a共x ជ ,t兲e ជ 储 .共14兲 Plugging this expression into the equation for the energy, we have the final result ⳵ E共x ជ ,t兲 ⳵ t⫽D⬜ⵜ⬜ 2 ␳ a共x ជ ,t兲⫺2␭ ⳵ 储 ␳ a共x ជ ,t兲 ⫺ ⑀ ␳ a共x ជ ,t兲⫹h共x ជ ,t兲⫹ ␩ E共x ជ ,t兲,共15兲 where the symbol ⳵ 储 stands for the partial derivative ⳵ / ⳵ x 储 . This is the general conservation equation for any directed sandpile model. It is worth remarking at this point that the energy field is a static field, in the sense that energy diffuses only if active sites are present in the system. This is intuitively understood in sandpile models, where energy 共sand兲 grains diffuse only from toppling sites. To analyze the consequences of Eq. 共15兲, it proves useful to define the susceptibility ␹ (x ជ ,t)关28兴: ␹ 共x ជ ⫺x ជ ⬘,t⫺t⬘兲⫽ 冓 ␦ ␳ a共x ជ ,t兲 ␦ h共x ជ ⬘,t⬘兲 冔 ␩ ,共16兲 where the symbol 具典 ␩ denotes an average over the noise distribution. By definition, the susceptibility measures the average increase in the number of active sites due to an impulsive perturbation, that is, to the addition of a single energy grain. Since we measure the size of the avalanches by the total number of topplings, the average avalanche size is given by 具 s 典 ⫽ 冕 ddxdt ␹ 共x ជ ,t兲.共17兲 Taking the functional derivative of Eq. 共15兲and averaging over time and noise, we obtain, in the limit t→⬁, in which the sandpile is in a stationary state with constant average energy, the following equation for the static susceptibility: D⬜ⵜ⬜ 2 ␹ 共x ជ 兲⫺2␭ ⳵ 储 ␹ 共x ជ 兲⫺ ⑀ ␹ 共x ជ 兲⫽⫺ ␦ (d)共x ជ 兲.共18兲 This equation can be easily solved in Fourier space. Defining the transformation ␹ 共x 储 ,x ជ ⬜兲⫽1 共2 ␲ 兲d 冕 dd⫺1kdq ␹ 共q,k ជ 兲eik ជ •x ជ ⬜eiqx 储 共19兲 and substituting into Eq. 共18兲, we obtain the solution ␹ 共q,k ជ 兲⫽1 D⬜k2⫹2i␭q⫹ ⑀ ,共20兲 which yields the susceptibility in real space ␹ 共x 储 ,x ជ ⬜兲⫽1 共2 ␲ 兲d 冕 dd⫺1keik ជ •x ជ ⬜ 冕 ⫺⬁ ⬁dq eiqx 储 D⬜k2⫹2i␭q⫹ ⑀ . 共21兲 This integral yields the result, setting D⬜⫽1: ␹ 共x 储 ,x ជ ⬜兲⫽1 2␭ 冉 ␭ 2 ␲ 冊 (d⫺1)/2 x 储 (1⫺d)/2e⫺x 储 ⑀ /2␭e⫺␭x⬜ 2/2x 储 . 共22兲 Equation 共22兲can be conveniently rewritten into the scaling form ␹ 共x 储 ,x ជ ⬜兲⫽x 储 (1⫺d)/2⌫ 冉 x 储 ␰ 储 ,x⬜ ␰ ⬜ 冊 ,共23兲 6202 PRE 62ROMUALDO PASTOR-SATORRAS AND ALESSANDRO VESPIGNANI where ⌫is a cutoff function that decreases exponentially in both its arguments. Comparing this last expression with Eq. 共22兲, we can identify the parallel and transversal correlation lengths ␰ 储 ⬃ ⑀ ⫺1, ␰ ⬜⬃ ⑀ ⫺1/2.共24兲 In more general terms, if we define the exponents ␯ 储 ⬘and ␯ ⬜ ⬘ by Eqs. 共11兲, then we have for directed sandpiles ␯ 储 ⬘⫽1 and ␯ ⬜ ⬘⫽1/2. From these last expressions, we can read off a first exact result for directed sandpiles: the avalanches produced in those models are elongated, with characteristic length in the parallel and transversal directions related by an affinity exponent ␨ ⫽1/2. It is very important to stress that these results are independent of the particular model considered and of the dimensionality dof the system, dictated only by the energy balance in the stationary state. We can use the result 共24兲to relate the effective bulk dissipation with the system size in a model with open boundary conditions. To sustain a steady state with constant average energy, avalanches must reach the bottom boundary in order to be able to dissipate. This means that the characteristic length of the avalanches in the parallel direction must be proportional to the system size ␰ 储 ⬃L. We have therefore that in boundary dissipation models we can define an effective dissipation rate ⑀ that is related with the system size by ⑀ ⬃L⫺1.共25兲 From this relations we easily find that ⌬s⫽Dand ⌬t⫽z. These identities are recovered in numerical simulations 共see Tables I and II兲. Finally, from Eq. 共20兲, we can recover the well-known result linking the system size and the average avalanche size, 具 s 典 ⫽ ␹ (q⫽0,k ជ ⫽0)⬃ ⑀ ⫺1⬃L关10,15,23兴. B. Anisotropic sandpiles Having completed the analysis of directed sandpiles, we turn our attention to the more complex case of anisotropic sandpiles. In this kind of model, the transport of energy is not strictly directed in the parallel direction, but is simply stronger in the direction ⫹x 储 than in the opposite direction ⫺x 储 . The presence of backwards flow allows the possibility of diffusion in the preferred direction, and thus the equation for the conservation of energy becomes in this case ⳵ E共x ជ ,t兲 ⳵ t⫽D⬜ⵜ⬜ 2 ␳ a共x ជ ,t兲⫹D 储 ⳵ 储 2 ␳ a共x ជ ,t兲⫺2␭ ⳵ 储 ␳ a共x ជ ,t兲 ⫺ ⑀ ␳ a共x ជ ,t兲⫹h共x ជ ,t兲⫹ ␩ E共x ជ ,t兲.共26兲 From Eq. 共26兲, we can obtain the corresponding equation for the susceptibility. The solution in Fourier space is readily found to be ␹ 共q,k ជ 兲⫽1 D⬜k2⫹D 储 q2⫹2i␭q⫹ ⑀ .共27兲 Upon integration over k ជ and q, one obtains the expression in real space ␹ 共x 储 ,x ជ ⬜兲⫽1 共2 ␲ 兲d 冕 dd⫺1keik ជ •x ជ ⬜ ⫻ 冕 ⫺⬁ ⬁dq eiqx 储 D⬜k2⫹D 储 q2⫹2i␭q⫹ ⑀ .共28兲 This last integral can be performed analytically in d⫽1 and 2. For d⬎2, even though we do not have a closed expression, we can obtain the leading scaling behavior. To simplify the calculations, we set, without lack of generality, D⬜⫽D 储 ⫽1. The integration in qis done by the method of the residues. The integration of the k ជ angular part 关38兴yields ␹ 共x 储 ,x ជ ⬜兲⫽1 2 冉 ␥ 2 ␲ 冊 ␯ ⫹1 x⬜ ⫺ ␯ 冕 0 ⬁dzz ␯ ⫹1 ⫻J ␯ 共 ␥ x⬜z兲e⫺x 储 ( ␥ 冑 1⫹z2⫺␭) 共1⫹z2兲1/2 .共29兲 Here, J ␯ (z) is the first kind Bessel function of order ␯ , and we have defined the constants ␯ ⫽(d⫺3)/2 and ␥ ⫽(␭2 ⫹ ⑀ )1/2. We are interested in the behavior of this integral for large distances, that is, in the limit x 储 Ⰷx⬜Ⰷ1. In this limit, the weight of the integral is given by the region of small z, since the exponential suppresses large values. We can then approximate the integral in the interval 0⬍z⬍1 and perform a Taylor expansion of the square root in the exponential and the denominator. In the denominator, we readily have (1 ⫹z2)1/2⯝1. The term in the exponential, however, contains a constant term, and must be therefore expanded up to second order: ⫺x 储 共 ␥ 冑 1⫹z2⫺␭兲⯝⫺x 储 共 ␥ 关1⫹z2/2兴⫺␭兲 ⫽⫺x 储 共 ␥ ⫺␭兲⫺x 储 ␥ z2/2. 共30兲 In the limit ⑀ →0, we have ␥ ⯝␭, and the constant ␥ ⫺␭can be expanded to give ␥ ⫺␭⫽共␭2⫹ ⑀ 兲1/2⫺␭⯝␭ 冉 1⫹ ⑀ 2␭2 冊 ⫺␭⫽ ⑀ 2␭.共31兲 Substituting these approximations into Eq. 共29兲, we are led to the expression ␹ 共x 储 ,x ជ ⬜兲⯝1 2␭ 1 共2 ␲ 兲 ␯ ⫹1x⬜ (1⫺d)e⫺x 储 ⑀ /2␭ ⫻ 冕 0 ⬁dy y ␯ ⫹1J ␯ 共y兲e⫺(x 储 /2␭x⬜ 2)y2,共32兲 where we have performed the change of variables y⫽ ␥ x⬜z and extended again the upper limit of the integral to infinity 共which is allowed given its exponential convergence兲. The integral in Eq. 共32兲yields 关38兴 ␹ 共x 储 ,x ជ ⬜兲⯝1 2␭ 冉 ␭ 2 ␲ 冊 ␯ ⫹1 x 储 (1⫺d)/2e⫺x 储 ⑀ /2␭e⫺␭x⬜ 2/2x 储 , 共33兲 PRE 62 6203CRITICAL BEHAVIOR AND CONSERVATION IN . . .