# Multifractal Scaling in the Bak-Tang-Wiesenfeld Sandpile and Edge Events

## Abstract

An analysis of moments and spectra shows that, while the distribution of avalanche areas obeys finite size scaling, that of toppling numbers is universally characterized by a full, nonlinear multifractal spectrum. Rare, large avalanches dissipating at the border influence the statistics very sensibly. Only once they are excluded from the sample, the conditional toppling distribution for given area simplifies enough to show also a well defined, multifractal scaling. The resulting picture brings to light unsuspected, novel physics in the model.

PACS numbers: 64.60.Lx, 64.60.Ak, 05.40.+j, 05.65.+b

Finite size scaling (FSS) is a widely adopted framework for the description of finite, large systems near criticality. In the last decade, after the work of Bak et al., much attention has been devoted to a class of models in which criticality is spontaneously generated by the dynamics itself, without the necessity of tuning parameters. This self–organized criticality (SOC) has been advocated as a paradigm for a wide range of phenomena, from earthquakes to interface depinning, economics and biological evolution . The prototype model of SOC is the two dimensional (2D) Bak, Tang, and Wiesenfeld sandpile (BTW) , which represents a system driven by a slow external influx, dissipated at the borders through a local, nonlinear mechanism.

In spite of its apparent simplicity and relative analytical tractability, the 2D BTW resisted, so far, all theoretical attempts to fully and exactly characterize its scaling. These attempts were essentially all based on the FSS ansatz. Numerical approaches, also based on FSS, led to rather scattered and sometimes contradictory numerical results , which hardly concile with existing theoretical conjectures. Thus, with its intriguing intractability, BTW scaling remains a formidable challenge for nonequilibrium statistical mechanics and it is very important to check if FSS works in this context.

In this Letter we apply a new strategy of data collection and interpretation, in order to determine to what extent the FSS ansatz can be applied, or rather has to be modified, for a correct description of the 2D BTW. Our results are striking and largely unexpected: while compelling evidence is obtained that the probability distribution functions (pdf) of some quantities obey FSS, for other magnitudes, whose fractal dimensions can widely fluctuate within the the nonlinear dynamics, this is definitely not the case. Following our protocol of analysis, we demonstrate that the well known difficulties in the description of the BTW are due to unexpected, very nonstandard features of its dynamical behavior. In the BTW, relations between different key quantities do not reduce to standard power laws, as in FSS, and are substantially influenced by the infrared cutoff given by the size of the system. The peculiar fluctuations characterizing intermittent dissipation at the cutoff scale, provide a dynamical mechanism for unusual deviations from finite size, and even multifractal scaling.

We consider the BTW on a square lattice box $`\mathrm{\Lambda }`$; to each site $`i`$ we associate an integer height $`z_i>0`$, the number of grains. When $`z_i`$ exceeds a threshold $`z_c=4`$, site $`i`$ topples: $`z_iz_i4`$, while for the nearest neighbors $`j`$ of $`i`$, $`z_jz_j+1`$. At the boundary less than $`4`$ neighbor sites are upgraded with consequent grain dissipation. Further instabilities can be created by the first toppling. An avalanche is the set of the $`s`$ topplings necessary to reach a stable system configuration after addition of one grain ( $`z_kz_k+1`$ at some randomly chosen $`k\mathrm{\Lambda }`$), $`a`$ is the number of lattice sites toppling at least once during the avalanche. A sequence of avalanches is created by successive random additions. After sufficiently many grains, thanks to dissipation at the borders, the sandpile reaches a steady state. We analyzed up to $`10^8`$ avalanches in this state for $`L=128,256,512`$ and $`1024`$.

The pdf for $`a`$ and $`s`$ do not reveal characteristic scales intermediate between lattice spacing and linear pile size $`L`$. FSS postulates for them the forms

$`P_s(s,L)`$ $`=`$ $`s^{\tau _s}F_s\left(s/L^{D_s}\right)`$ (1)   
$`P_a(a,L)`$ $`=`$ $`s^{\tau _a}F_a\left(a/L^{D_a}\right),`$ (2)   
and usually assumes that different quantities characterizing an avalanche are simply related by power laws, determined by sharply peaked conditional pdf. For example, one expects that, given $`a`$, $`sa^\gamma `$ , with $`D_s=\gamma D_a`$, and $`D_s`$ ($`D_a`$) representing the fractal dimension of $`s`$ ($`a`$). Without sticking to the FSS form (1), scaling of $`P_s`$ is most generally described by the multifractal spectrum:

$$f\left(\alpha \right)=\frac{\mathrm{log}(_{L^\alpha }^{\mathrm{}}P_s(s,L)𝑑s)}{\mathrm{log}(L)}.$$
(3)   
$`f`$ is the Legendre transform of the moment scaling function $`\sigma \left(q\right)`$ defined by:

$$s^q_L=P_s(s,L)s^q𝑑sL^{\sigma \left(q\right)},$$
(4)   
i.e. $`\sigma \left(q\right)=sup_\alpha [q\alpha +f(\alpha )]`$. Analogous definitions apply to the spectrum $`g(\beta )`$ ($`\beta =\mathrm{log}a/\mathrm{log}L`$) and the moment exponent $`\rho (q)`$ of $`P_a`$. If Eqs. (1) hold, $`f\left(\alpha \right)=(\tau _s1)\alpha `$ for $`0<\alpha <D_s`$ and $`f=\mathrm{}`$ for $`\alpha >D_s`$. Consistently $`\sigma \left(q\right)=D_s(q\tau _s+1)`$ for $`q>\tau _s1`$ and $`\sigma \left(q\right)=0`$ for $`q<\tau _s1`$. Corresponding expressions for $`a`$-quantities hold if FSS is valid. So, within FSS, both $`f`$ and $`\sigma `$ are piece-wise linear functions of their arguments.

A reliable way to establish if Eqs.(1) hold, is by checking the above linearity of $`\sigma `$ and $`\rho `$ in significant ranges of $`q`$. Being extrapolated for $`L\mathrm{}`$ from finite $`L`$ moments, $`\rho `$ and $`\sigma `$ provide a very asymptotic characterization of the pdf’s. While for $`P_a`$ a constant gap $`\mathrm{\Delta }\rho (q)=\rho (q+1)\rho (q)2.02\pm 0.03D_a`$ establishes already for $`q=1`$ (Table I), for $`P_s`$, $`\sigma ^{}(1)2.5`$, while $`\mathrm{\Delta }\sigma `$ steadily increases from $`\mathrm{\Delta }\sigma (1)2.70`$ to $`\mathrm{\Delta }\sigma (8)2.92`$. $`\mathrm{\Delta }\sigma (1)2.7`$, was also found in Ref. , based on large $`L`$ data. Thus, unlike for $`P_a`$, and in violation of FSS, for $`P_s`$ there is clearly no constant gap in the range $`q1`$. The gap tends to rise to $`3.0`$ for increasing $`q`$. Therefore, we should also expect $`f(\alpha )>\mathrm{}`$ as long as $`\alpha 3.0`$.

Fig.1 reports $`g`$ and $`f`$ as obtained by the data collapse technique in Ref. . The linear form $`g\left(\beta \right)=(\tau _a1)\beta `$ is well verified for $`\beta 1.5`$ with an estimated slope $`0.19\pm 0.01`$. In the region $`\beta >1.5`$, the collapse gets worse. One expects $`D_a=2`$ and, in fact, the poor $`g`$-collapse for $`\beta 2`$ is consistent with the infinite discontinuity of a FSS spectrum with $`D_a=2`$: curves for various $`L`$ smooth out to different degrees such discontinuity, and underestimate $`g`$ for $`\beta 2`$. Assuming a linear $`g`$ in the whole domain $`0\beta 2`$ , and $`\tau _a=1/5`$ exactly, as suggested by the estimated initial slope, $`a_L`$ should scale with $`\rho (1)=sup_\beta \left[g(\beta )+\beta \right]=2/5+2=1.6`$, in nice agreement with our determination $`\rho (1)=1.59\pm 0.02`$. Evidence of FSS for $`P_a`$ comes also from the fact that a standard FSS collapse ($`\tau _a=6/5`$; $`D_a=2`$) works very well (Fig. (1), inset). For $`\alpha 2`$ the collapsed $`f`$ is very close to linear and overlaps with the expected $`g`$ ($`f\left(2\right)0.39`$). Thus, an acceptable FSS form of $`P_s`$ should assume $`\tau _s=\tau _a`$, in order to be consistent with the well collapsed, initial part of the plots. Within such assumption, $`D_s=2.5`$ would be imposed by the exact result $`\sigma (1)=2`$ (we find $`\sigma \left(1\right)=1.99\pm 0.02`$). Indeed, $`sup_\alpha [f(\alpha )+\alpha ]=2`$ in such case, the $`sup`$ being attined at $`\alpha =D_s=2.5`$. The hypothetical linear spectrum should therefore have support in $`0\alpha 2.5`$, and satisfy $`f\left(5/2\right)=1/2`$. In Ref. , a relatively limited analysis of $`\sigma `$ suggested $`\mathrm{\Delta }\sigma (q)2.5`$ for $`q1`$. Such constant gap would leave room to a FSS approximation for $`P_s`$, with $`\tau _s6/5`$ and $`D_s2.5`$, and a linear $`\sigma `$ deviating from the measured one possibly only at very low $`q`$’s. According to Eq. (1), for such effective $`P_s`$ one would also have $`D_s(2\tau _s)2=\sigma (1)`$. However, this does not agree with the plots which show that $`f>\mathrm{}`$ at least up to $`\alpha 3.0`$. Indeed, curves for various $`L`$ collapse rather well for $`\alpha 3.0`$, and clearly suggest a support of $`f`$ asymptotically bounded by $`\alpha 3.0`$, consistent with the trend of $`\mathrm{\Delta }\sigma `$ (Table I). This bound follows also from the leftward trend of the curves for increasing $`L`$ in the region $`\alpha 3`$, where collapse gets worse. An even approximate $`\tau _s`$ exponent, so extensively discussed in the last decade, can not be simultaneously consistent with the initial slope $`1/5`$ of $`f`$, $`\mathrm{\Delta }\sigma (\mathrm{})3`$ and $`\sigma (1)=2`$. Full consistency, in particular with $`\sigma (1)=2`$, can be recovered by assuming that $`f`$ is indeed linear, with the same slope $`1/5`$ as $`g`$, up to $`\alpha =2.5`$ ($`f(2.5)=1/2`$), but has a nonlinear continuous drop in the range $`2.5<\alpha <3.0`$. The slight underestimation by the plot for $`\alpha 2.5`$ ( $`f(5/2)0.57`$) should again be imputed to roundoff and slower $`L`$-convergence in correspondence to the major bending of $`f`$ .

The striking difference described above between $`P_s`$ and $`P_a`$ suggests that the conditional pdf, $`C`$, such that

$$P_s(s,L)=𝑑aC(s|a,L)P_a(a,L),$$
(5)   
should have unusual structure. Within FSS one would assume $`C\delta (sa^\gamma )`$ for $`a<L^2`$ and $`L\mathrm{}`$. Such $`C`$ would lead to $`\gamma =(\tau _s1)/(\tau _a1)`$. Thus, $`\tau _s=\tau _a`$ would imply $`\gamma =1`$, while $`D_s=2.5`$ would give $`\gamma =D_s/D_a=1.25`$, already contradicting FSS. In fact $`C`$ is a complex, broad pdf. Its properties elucidate and confirm our conclusions for $`P_a`$ and $`P_s`$. Fig. 2 reports values of $`\alpha =\mathrm{log}s/\mathrm{log}L`$ vs. $`\beta =\mathrm{log}a/\mathrm{log}L`$ for $`L=1024`$ avalanches. Ratios $`\gamma =\alpha /\beta `$ range between $`1`$ ( $`sa`$) and $`1.25`$, and their spread is not modified appreciably by sampling data for progressively larger $`L`$’s. For $`\beta 2`$ also some $`\alpha /\beta >1.25`$ are found, which are too rare to be displayed in Fig. 2. If a relation $`sa^\gamma `$ would hold, data should coalesce into a straight line with slope $`\gamma `$. To the contary, points are quite spread and form an open angle, rather than a narrow strip, as is instead the case, e.g., for the radius of gyration, $`r`$, of the surface covered by the avalanche (Fig. 2).

$`C`$-moments can be assumed to scale as $`s^q_a=𝑑sC(s|a,L)s^qa^{\kappa (\beta ,q)}=L^{\beta \kappa (\beta ,q)}`$. One would hope exponents like $`\kappa `$ to be well defined, i.e. independent of $`\beta `$, for $`L\mathrm{}`$. Furthermore, with a FSS $`C\delta (sa^\gamma )`$, one should find $`\beta \kappa (\beta ,q)/q=\beta \gamma `$ independent of $`q`$. In Fig.3 we plot $`\mathrm{log}(s^q_a^{1/q})/\mathrm{log}L`$ versus $`\beta `$ for various $`q`$ and $`L=1024`$. The curves, which correspond to $`\beta \kappa (\beta ,q)/q`$ asymptotically, do not overlap. Moreover, for each $`q`$ the plots have pronounced curvature, especially for $`\beta 1.7`$. This curvature does not decrease appreciably, or increase signalling crossovers, if one considers progressively larger $`L`$’s. Thus, in spite of the sensible spread of the various curves in Fig.3 (which also persists for increasing $`L`$), the complex scaling of $`C`$-moments can not even be classified as multifractal. Indeed, in that case one should still require $`\beta \kappa /q`$ to be linear in $`\beta `$. On the other hand, the $`\beta `$-nonlinearity described above is essential in view of our conclusions on $`g`$. Since $`s^q_L=𝑑aP_a(a,L)s^q_a`$, we must have $`\sigma (q)=sup_\beta [g(\beta )+\beta \kappa (q,\beta )]`$. With $`g=\beta /5`$, as argued above, the $`sup`$ for $`q=1`$ should fall at $`\beta =2`$ (Fig. 3), and be equal to $`2/5+2\kappa (1,2)`$. Thus, one must find $`\kappa (1,2)=1.2`$. If $`\kappa (1,\beta )`$ would remain constant with its initial value $`1.04\pm 0.02`$ in the whole $`\beta `$-range, we could certainly not find $`\sigma (1)=2`$. Our determinations consistently extrapolate to $`\kappa (1,2)1.20`$ (Fig. 3).

The $`\beta `$-dependence of $`\kappa `$ is crucial for the consistency of the overall scaling properties and is determined by large dissipating avalanches. These are intermittent, edge events, which guarantee grain conservation at stationarity. Remarkably, if only nondissipative avalanches are sampled, one obtains a different conditional pdf, $`C_0`$, whose $`\kappa _0`$ is now independent of $`\beta `$ (Fig.3). Thus, for $`C_0`$ it makes sense to discuss a multifractal spectrum

$$\beta h(\gamma )=\frac{\mathrm{log}(_{L^{\beta \gamma }}^{\mathrm{}}𝑑sC_0(s|a,L))}{\mathrm{log}(L)}.$$
(6)   
As a function of $`\gamma `$, $`h`$ measures the density of points ($`\alpha =\beta \gamma ,\beta `$) along a narrow vertical strip at fixed $`\beta `$ in the plot correspoding to that in Fig.3. Fig.4 reports a data collapse of $`h`$-curves at $`\beta =1.6`$. Even if the collapse suffers of poor $`L`$-asymptoticity, all curves cross nicely at a point corresponding to $`\gamma 1.25`$, where their trend for increasing $`L`$ inverts itself. Collapses at different $`\beta `$’s give very similar results. Thus, $`\gamma 1.25`$ qualifies as a maximum exponent for nondissipating, bulk avalanches. An independent, accurate determination based on rank ordering statistics gave $`\gamma =1.28\pm 0.03`$.

$`\gamma =5/4`$ was conjectured as a unique exponent for all avalanches in Ref.. Different values, sometimes larger than $`1.25`$, were conjectured in the recent literature on the basis of FSS and numerical results. Figs. 3 and 4 show that the determination of $`\gamma `$ is not even a well defined task, and can be discussed, within a multifractal framework, only once edge avalanches are eliminated. Indeed, for nondissipative bulk avalanches, $`\gamma `$ has a $`\beta `$-independent, broad, nonlinear spectrum, $`h`$, with $`1`$ as most probable, and $`1.25`$ as extreme, most rare values.

The above analysis identifies the inadequacies of the scaling theory in Ref.. Crucial to that approach was the assumption of a relation, $`na^{(2\gamma 2)/2}=a^{1/4}`$, satisfied by $`n`$, the maximum number of topplings in an avalanche (realized at the injection point), and $`a`$, seen as the area of the first in a succession of waves . We found that also the $`n`$-pdf at given $`a`$ is broad, and $`1/4`$ is only the approximate maximum value of the exponent realized by bulk sandpile dynamics.

The constant gap $`\mathrm{\Delta }\sigma 2.5`$ at high $`q`$ postulated in Ref.  would have been compatible with a narrow, or even point-like spectrum for dissipating avalanches ($`\alpha 2.5`$). To the contrary, it turns out that also dissipating avalanches are multifractal ($`2.0\alpha 3.0`$) consistent with our more correct estimate of $`\mathrm{\Delta }\sigma `$ at high $`q`$. Moreover, a frequency $`L^\zeta `$ with $`\zeta =\frac{1}{2}`$ is precisely associated to dissipating avalanches with $`\alpha 2.50`$.

In summary, we showed how the long standing puzzle of 2D BTW scaling finds a solution in a strong violation of FSS by both $`P_s`$ and $`C`$. While $`P_a`$ obeys FSS, with $`D_a=2`$ and $`\tau _a=6/5`$ as most plausible exponents, nonlinear multifractal spectra are needed to characterize $`P_s`$ and $`C_0`$. Our results throw light on the intriguing difficulty of this model. In spite of the several exactly known steady state properties, the belief that the 2D BTW should be “easily” solvable appears unjustified. The unusual scaling pattern discovered here, which is not found in more simple systems, like directed sandpiles , constitutes genuine novel physics and enhances the paradigmatic role of the BTW for statistical mechanics out of equilibrium.

The most striking feature of this pattern is the bending of the curves in Fig.3, determined by the effect of intermittent, edge avalanches, and essential in order to fulfill the exact constraint of local, Laplacian conservation at stationarity ($`\sigma (1)=2`$). In fact the very presence of edge avalanches allows the existence of a peculiar, bulk multifractal scaling of $`C_0`$. The first moment of the $`s`$-pdf of nondissipative avalanches scales $`L^{1.8}`$, rather than $`L^2`$, as conservation imposes in the case of $`P_s`$. The role played here by intermittent edge avalanches has striking analogies with the intermittent phenomena in fully developed turbulence, where they cause deviation from pure Kolmogorov scaling. Results of the present analysis have been used most recently to elucidate quantitative connections between BTW scaling and the notoriously difficult problem of turbulence in 3D .

We acknowledge partial support from the European Network Contract No. ERBFMRXCT980183. We are grateful to D. Dhar for useful criticism.