Roughness of Sandpile Surfaces

We study the surface roughness of prototype models displaying self-organized criticality (SOC) and their noncritical variants in one dimension. For SOC systems, we find that two seemingly equivalent definitions of surface roughness yields different asymptotic scaling exponents. Using approximate analytical arguments and extensive numerical studies we conclude that this ambiguity is due to the special scaling properties of the nonlinear steady state surface. We also find that there is no such ambiguity for non-SOC models, although there may be intermediate crossovers to different roughness values. Such crossovers need to be distinguished from the true asymptotic behaviour, as in the case of a noncritical disordered sandpile model studied in [10].

Since the original proposal by Bak, Tang and Wisenfeld [1] there has been a large body of work directed towards understanding ubiquity of scale invariance in externally driven open dissipative systems using the concept of self-organized criticality (SOC). The sandpile is a prototype model system which has been extensively used as a paradigm of SOC [2,3]. The principal aim has been to elucidate how slowly driven dissipative systems with fast relaxation mechanisms display long-tailed distributions of activity sizes. Recently, there have been resurge of interest in sandpile models in order to understand SOC in connection with better understood scale invariance in other nonequilibrium phenomena such as absorbing state phase transitions [4] and driven interfaces [5].
In [6], Krug, Socolar and Grinstein (KSG) studied the surface fluctuations in a prototype model of SOC, the limited local sandpile (LLS) [3] and its variations in one spatial dimension and concluded that the interfacial fluctuations, although nontrivial, are evidently unconnected to the criticality of the system. In fact, they argued that the evolution of height fluctuationsh(x, t) can be described by an extension of the Kardar-Parisi-Zhang (KPZ) [7,8]

equation for an anchored interface
where η(x, t) is a Gaussian white noise. The linear term c∂ xh is more relevant than the diffusion term D∂ 2 xh and the KPZ nonlinearity λ(∂ xh ) 2 (in RG sense) and, since the interface is anchored, can not be eliminated by a Galilean shift. This term is responsible for transporting fluctuations up the pile thus relating the spatial roughness of the anchored interface to the dynamical roughening of a moving KPZ interface in d = 1, α LLS = β KP Z(1d) = 1/3. This argument can be easily extended to higher dimensions. In this paper we study interface roughness in LLS and related models and find that effect of the criticality of the system is somewhat subtle and is due to the nontrivial nature of the steady state surface. The value of surface roughness depends on whether it is measured about this nontrivial steady state profile or about the mean instantaneous surface, the latter being the correct choice for moving interfaces. The assertion made by KSG is valid only ifh in (1) is defined as the fluctuation around the nonlinear steady state surface. We also study modifications of the LLS which have nonlinear steady state profiles but which do not display SOC and find that there is no such ambiguity in the asymptotic roughness exponent. Thus, we conjecture that the roughness exponent is uniquely defined for systems not displaying SOC even though they may possess nontrivial steady state surfaces [9]. For SOC systems the ambiguity exists and one needs to define the roughness appropriately. Our results have important bearings on studies such as the recent one by Barker and Mehta [10] who observed anomalously large roughness exponent in a sandpile model with structural disorder.
The limited local sandpile model in d = 1 is defined as follows (see Fig. 1a): on a one dimensional lattice of sites i = 0, 1, · · · , L we define an integer height variable h i . One grain of sand added at a randomly chosen site increases the height at that site by one: The configuration is stable if the local slopes satisfy where we chose z c = 2 (our results are essentially unchanged for z c ≥ 2). An instability occurs when by addition of grains at site i + 1, the local slope exceeds threshold z i > z c ; in this case z c grains are transferred from column i + 1 to column i. Subsequently, columns i and i + 2 may become unstable, setting off further topplings leading to an avalanche. The grains involved in an avalanche leave the system if they reach site i = 0 (h 0 (t) = 0, ∀t). A new grain is added after the system has attained a stable configuration. The final configuration reached is independent of the order in which the sites are updated in case more than one site become unstable. We count one unit of time (Monte Carlo step) for every L grains added.
The standard way of quantifying avalanches in these systems is to count the number m of grains that leave the pile after each deposition [3]. The criticality of the model is reflected in a broad probability distribution of avalanche sizes, P (m) (Fig. 1b), which has no other length scale except the system size L.
In [6], KSG also introduced and studied the inertial limited local sandpile (ILLS) which mimics the effect of inertia of the falling grains in a real sand pile. The instability condition setting off an avalanche is same as that in the LLS but the condition of stopping is changed. A cluster of grains when first destabilized is assigned a momentum p = 0 and each time the front of the cluster reaches an unstable site (z i > z c ) it gathers momentum, p → p + 1. If it comes across a stable site, p decreases by an amount r: p → p − r. The cluster continues moving down as long as p > 0 and leaves the pile if it reaches the site i = 0. Clearly, r = ∞ corresponds to the LLS. It was noted by KSG that ILLS is not critical for any finite value of r which is reflected in the corresponding avalanche size distribution (Fig. 1b).
To study the interfacial fluctuations in these models, we start with an initially flat pile (h i (0) = 0, ∀i) and add grains till the pile reaches a steady state with the mean surface making a critical angle ψ (tan ψ L→∞ = 3/2) with the horizontal. The closed boundary condition at i = L ensures that the pile has only one surface with the critical slope [9]. The width of the interface, in the steady state, can be measured in two ways where In both definitions (2) above the ensemble average · · · is identical to the time average in the steady state. For large enough system sizes the scaling hypothesis is expected to be valid and the roughness exponent α is defined through In Fig. 2, we plot the two widths W 1 and W 2 for a range of system sizes L = 2 n with 6 ≤ n ≤ 14 (that is 64 ≤ L ≤ 16384). It is seen that the two widths have different asymptotic behaviours and give two different values of α. The roughness α ≃ 0.33 computed from W 1 (and W 2 up to L < ∼ L c ) is in accordance with the predictions of Eq.(1). The asymptotic roughness computed from W 2 , beyond the crossover system size L c ∼ 1024, α = 0.65 ± 0.02 is larger than that for an interface generated by a simple random walk. This is surprising since, prima facie, the definitions in (2) are expected to be equivalent as far as asymptotic scaling is concerned ( [11]). In fact, for the ILLS, we have computed W 1 and W 2 for r = 20 (Fig. 2, inset) and it is seen that they indeed have the same asymptotic behaviour and hence, a unique value of α ≃ 0.33.
We now show that the crossover in W 2 in the LLS can be traced to the fact that the time averaged steady state profile is not linear, i.e., h i (t) =s i and this makes the two definitions nonequivalent. Although, this nonlinearity of h i (t) is a consequence of the boundary effects and is present for the ILLS as well, we demonstrate that the presence of SOC in LLS results in special scaling properties which is responsible for the observed differences between asymptotic scalings of W 1 and W 2 . It can be shown that W 1 and W 2 are related by where W S (L) is the root-mean-square (rms) wandering of the steady state interface profile h i , withs = s(t) being the mean slope of the steady state profile satisfying L i=0 [ h i −s · i] = 0. The second term on the left of (4), W 2 c (L) = [s(t)−s] 2 (L+1)(2L+1)/6, is the contribution to interface width due to instantaneous slope fluctuations and is the analogue of the center of mass fluctuations in the case of moving interfaces [8]. In fact, in the case of moving interfaces W c dominates for large length scales and hence W 2 is taken as the width. As it turns out, for the anchored interfaces we deal with here, W c is negligible compared to the other terms in (4) as L → ∞. From Fig. 2, it is evident that the crossover in W 2 (L) for the LLS is due to the contribution from W S . In the following we show that, for the LLS, the special form of the nonlinear steady state profile h i (t) , which results in W S ∼ L 0.66 , is a consequence of the singular diffusion associated with the SOC state.
In [12], using a simple model, Carlson et. al. demonstrated that the SOC state of the system is associated with the vanishing of the density of troughs concommittant with the divergence of the corresponding diffusion coefficient. The troughs are defined as the sites for which z i ≤ 0 so that, in the LLS, an avalanche necessarily stops at a trough (or leave the system at i = 0). On a coarse grained level we define a set of densities {ρ n (x, t); n = −∞, · · · , −1, 0, 1, 2} where ρ n (x, t) denotes the local density of sites with z i = n. It follows that the coarse grained local slope z(x, t), which is locally conserved by the dynamics, may be expressed in terms of the ρ n 's as z(x, t) = 2 n=−∞ nρ n (x, t). The open boundary condition at i = 0 implies h(0, t) = 0 and the closed boundary condition at i = L is modeled by setting z(L, t) = 0. The trough density ρ(x, t) ≡ 0 n=−∞ ρ n (x, t) is not strictly conserved. In the '01' model considered in [12] the slope has only two values z i = 0 (trough),1 and hence both z as well as ρ are strictly conserved and are related simply as z = 1 − ρ. Although ρ is not strictly conserved for LLS and ILLS we still approximate its dynamical evolution by the continuity equation, ∂ t ρ(x, t) + ∂ x J(ρ(x, t)) = 0. Phenomenologically, the trough current J consists of three parts: (i) current due to addition of grains J 0 , (ii) avalanche current J A , and (iii) a microscopic noise term η(x, t). In analogy with driven diffusive systems [6] the phenomenological form of the avalanche current is written as J A (ρ) = aρ + bρ 2 + · · · − D(ρ)∂ x ρ. As in the 01 model, the critical state of the system (L → ∞, r = ∞) is associated with ρ, ∂ x ρ → 0 and hence in order to balance the finite input flux J 0 , the diffusion constant D(ρ → 0) has to diverge appropriately [12]. Without loss of generality the leading divergence in D is taken to be a simple pole and the evolution of ρ(x, t) is thus [13], Here A(ρ) is taken to be a smoothly varying function in the relevant interval 0 ≤ ρ ≤ 1.
In the following we find an approximate numerical relation between ρ(x) to z(x) which would enable us to find the steady state interface profile h 0 (x) = x 0 z(x)dx form (7). In Fig. 3 (inset), we plot the average densities ρ n ≡ L −1 L 0 dx ρ n (x, t) in the steady state as functions of system size L. We note that, as L increases, while ρ 1 ≃ ρ 2 ≈ 1/2 remain finite, densities of all the trough sites vanish algebraically with L: ρ 0 , ρ −1 ∼ L −0.33 , ρ −2 ∼ L −0.6 , and ρ n 's with n ≤ −3 are negligible. Hence, in the limit of large L, we can approximate the total trough density as ρ ≃ ρ 0 + ρ −1 . From the normalization condition 2 n=−∞ ρ n (x, t) = 1 it follows that ρ 1 + ρ 2 = 1 − ρ and we find numerically that ρ 1 − ρ 2 ≃ 0.46ρ. Thus, for large systems, we may write z(x) ≃ 2ρ 2 + ρ 1 − ρ −1 ≃ 3/2 − κρ(x), where numerically κ ≃ 2 [14]. Thus, h 0 (x) is given approximately by In Fig. 3, using the approximate expression for ρ(x) from (7) with φ = 4, we compare h 0 (x) with that obtained numerically and notice that the agreement is rather good given the nature of approximations involved [15]. The L dependence of W S of the steady state profile can be computed from (8) and turns out to be For the LLS with θ ≃ 0.33, W S (L) ∼ L 0.66 dominates W 1 ∼ L 1/3 and W c ∼ L 0.18 in (4) and hence, for large L, W 2 (L) ∼ L 0.66 -in very good agreement with our numerical results (Fig. 2). It is interesting to note that W 1 and W 2 would have different asymptotic behaviours if θ < 2/3, i.e., if φ > 5/2. For the '01' model studied in [12] it is shown that φ = 3 exactly and thus W 2 (L) ∼ L 1/2 for the corresponding interface [16]. Our conjecture that W 1 and W 2 would have different asymptotic behaviour for SOC systems would be invalid if one can devise a model with φ < 5/2. Next, we present a general argument as to why, for noncritical models such as the ILLS, we do not expect any ambiguity in roughness exponent. The essential difference between a critical and a noncritical model is the presence of an additional length scale (apart from system size and the microscopic cutoff) which shows up in the avalanche size distribution (e.g. Fig. 1 for ILLS) as well as governs the decay of the boundary effects into the bulk. In the ILLS this length scale is related to the parameter r. To see this we first note that as a cluster moves down its momentum p makes a random walk [6]. If the mean density of troughs ρ < 1/r then most avalanches leave the system resulting in net drainage and if ρ > 1/r then most avalanches stop on the pile leading to net growth of the pile. Thus, in the steady state one has a finite density of troughs ρ = 1/r and thus the mean spacing between the troughs 1/ρ = r appears as the additional length scale. As already discussed above, the current J(ρ) now has a nonzero systematic part aρ + bρ 2 . . ., in addition to the finite diffusion term. Such a systematic current, e.g., a term such as aρ, results in the effects of the boundaries decaying exponentially inside the bulk. This is in contrast to the long ranged power law decay in the LLS which is reflected in the forms of ρ(x) and h 0 (x). Although the steady state surface is nonlinear (Fig. 3), its rms wandering W S is bounded and does not alter the asymptotic scaling of W 2 .
Lastly, we briefly point out possible pitfalls of using W 2 naively without properly accounting for the underlying steady state profile, even in systems which do not show SOC. Recently, Barker and Mehta [10] studied a disordered version of the LLS where disorder was introduced by allowing the added grains to have an aspect ratio different from unity: grains are now rectangles and are deposited on the sandpile with fixed probabilities of landing on their larger or smaller edges. Thus, the height h i of column i, no longer an integer, is the sum of the vertical edges of all the grains in that column. In addition to the threshold dynamics of LLS, dynamical reorientation of cluster of grains were allowed. They found that the larger sandpiles cease to display SOC which is reflected in the emergence of a preferred size of the large avalanches in the associated drop number distributions. Their numerical studies of W 2 showed that while for very small systems (L < ∼ 100) the roughness exponent α ≃ 0.34, it seems to crossover to a much larger value α ≃ 0.72 for larger system sizes, 100 < ∼ L < ∼ 400 (Fig. 4). Our preliminary numerical studies, using the same parameters as in [10], of yet larger systems (500 < ∼ L ≤ 2048) show that in fact the crossover seen in [10] is evidently transient and the asymptotic roughness exponent goes back to α ≃ 0.33 (see Fig. 4). Interface width of the disordered sandpile as a function of system size. The squares represent W2 and circles represent W = W 2 2 − W 2 S . The parameters chosen are the same as in [10]. The straight lines have slope 0.33 (dotted) and 0.72 (dashed). Inset: The steady state interface profile, after subtracting the mean linear profile.
In fact, the crossover reported in [10] is less noticeable if one looks at W = W 2 2 − W 2 S instead of W 2 (Fig. 4) [17]. This is in accordance with our conjecture since the disordered system is not critical, although the steady state profile is nonlinear (Fig. 4, inset), it does not affect the asymptotic roughness measured by W 2 . The transient crossover seen in W 2 here is similar to what is seen for the ILLS in Fig. 2, inset.
In summary, we have studied the surface roughness of a prototype model of SOC and its modifications in one dimension. We find that one needs to be careful in defining quantities such as the interface width since special form of the steady state shape of the surface in systems with SOC can result in different asymptotic behaviours of otherwise equivalent definitions. Although there is no such ambiguity for noncritical models, still there may be crossovers at intermediate length scales which should not be taken as the true asymptotic behaviour.
JGO and GT acknowledge financial support of Centro de Física do Porto, University of Porto where part of this work was carried out. JGO also acknowledges financial support of Universidade de Aveiro. JFFM was partially supported by project POCTI/1999/FIS/33141 and POCTI/2002/MAT/46176.