Directed Abelian sandpile with multiple downward neighbors

We study the directed Abelian sandpile model on a square lattice, with $K$ downward neighbors per site, $K>2$. The $K=3$ case is solved exactly, which extends the earlier known solution for the $K=2$ case. For $K>2$, the avalanche clusters can have holes and side-branches and are thus qualitatively different from the $K=2$ case where avalance clusters are compact. However, we find that the critical exponents for $K>2$ are identical with those for the $K=2$ case, and the large scale structure of the avalanches for $K>2$ tends to the $K=2$ case.


I. INTRODUCTION
The directed Abelian sandpile model is a simple variation of the sandpile model first introduced by Bak, Tang and Wiesenfeld [1]. It was the first model of selforganized criticality solved exactly [2], and its critical exponents can be determined in all dimensions d. The exponents take the classical values for d > 3, and for d = 3, there are logarithmic corrections to power-law behavior [3]. The model is related to other models of non-equilibrium statistical physics like the voter model, Scheideggar river network model, and Takayasu model of aggregating and diffusing particles with injection [4]. It has also found applications in more complex situations like modelling economic networks [5], growth of droplets of water in falling rain [9] and fracture of ice-sheets [7].
While the exact solution of the directed Abelian sandpile model on a square lattice with K = 2 downward neighbors per site is rather elementary, it depends crucially on the fact that avalanche clusters in this case have no holes and thus the problem reduces to that of two annihilating random walkers. If we consider sandpile models with K > 2, this property is no longer true, and it is not clear if the problem for K > 2 belong to the same universality class as K = 2. In fact, direct estimates of critical exponents from Monte Carlo simulations of the model with K > 2 show a persistent deviation from the exactly known K = 2 values [8].
The aim of this paper is to resolve this discrepancy. We provide an exact solution for the K = 3 case, and show that both K = 2 and K = 3 belong to the same the universality class. While this conclusion is not very surprising, the exact solution for K = 3 is of some in-terest, as we get the exact expression for the generating function of the mean-squared flux at a given depth from the top. This was not done in the earlier study of the K = 2 case, where only the critical exponents were deduced using scaling arguments. The K = 3 avalanches differ from K = 2 avalanches mostly near the surface, as can be seen from the pictures of typical avalanche clusters ( Fig. 1). We find that the presence of holes and side-branches in the avalanches for K = 3 is irrelevant in the sense of renormalisation group theory. Finally, we present the results of large scale Monte Carlo studies of this model for K = 2, 3, and 4. Our data is consistent with all of these being in the same universality class, and the observed deviations of exponents from the exact theoretical values in earlier simulations may be ascribed to the effect of significant corrections to scaling.

II. DEFINITION OF THE MODEL
We consider a directed sandpile model on a square lattice. The lattice is of size (height) L and width M . The sites are labelled X ≡ (x, t), where x ∈ {0, . . . , M − 1} and the vertical coordinate, thought of as time coordinate, is denoted by t ∈ {0, 1, . . . , L − 1}. The top row is t = 0 and the t-coordinate increases downwards (Fig.  2). We assume periodic boundary conditions in the xdirection, so that the x-coordinate is defined modulo M .
At each site X ≡ (x, t), there is a non-negative integer variable h(x, t), called the height of the pile. If h(x, t) < K the site is said to be stable. If h(x, t) ≥ K, the site is unstable and is said to topple. When K = 3,  on toppling at (x, t), it will send one particle each to its three downward neighbours (x − 1, t + 1), (x, t + 1) and (x + 1, t + 1).
Since the avalanche propagation depends only on the layers below the site where the particle is added, without loss of generality, we may assume that particles are added only in the top layer. Thus, we add particles at randomly selected sites in the row t = 0 and let the configuration relax by toppling. If a toppling occurs at any site in the bottom layer t = L − 1, three particles are lost from the system. The total number of topplings after adding a particle to the top row is called the avalanche size s and the set of sites where topplings occur, the corresponding avalanche cluster.
One can easily demonstrate that this is an Abelian[? ] model. There are K LM stable recurrent configurations and all occur with equal probability in the steady state, if particles are added everywhere [4]. However, if particles are added only in the top layer, then the system breaks into K L−1 disjoint sectors. In the thermodynamic limit of large lattice, i.e., L, M 1, the avalanche statistics in the steady state is independent of the sector.

III. CALCULATING THE 2-POINT CORRELATION FUNCTION
Consider the pile in the steady state and add a particle at the origin O ≡ (0, 0). Let η( X) denote the indicator variable that this causes a toppling at the site X ≡ (x, t). We define the two points correlation function G 2 (x, t) as the expected number of topplings at X, given that O topples after adding a particle at O, i.e., Clearly, G 2 (0, 0) = 1 as we are considering the expected no. of topplings under the condition that site O topples.
In the steady state, any given site X = (x, t) has equal probability of being with height 0, 1 or 2. Thus, if we focus on the activity in the layer above, the probability of topplings at (x, t) in a single avalanche equals to 1/3 times the expected number of upward neighbors that have toppled in the avalanche. Thus, G 2 (x, t) satisfies the equation With the boundary condition G 2 (x, t = 0) = δ x,0 , this equation determines the function G 2 (x, t) for all x and t.
If we define the characteristic functioñ it is easily seen that Substituting this into an inverse transform of Eq. (3), we get Note that G 2 (x, t) vanishes strictly for all |x| > t. We call the region |x| ≤ t the future light cone of the origin, see Fig. 1. For large t, G 2 (x, t) is well-approximated by a Gaussian of zero mean and variance 2t/3. We define Φ(t) as the random variable that measures the number of topplings in the layer t caused by driving at O: Then, it is easily seen that We would now like to calculate the mean square flux Φ 2 (t) . By definition: where Note that since each site topples at most once, G 3 (x 1 , x 2 |t) is in fact the probability that sites X 1 and X 2 both topple in a given avalanche.
Let us denote (x 1 + δ 1 , t − 1) and (x 2 + δ 2 , t − 1), δ 1 , δ 2 ∈ {−1, 0, 1} to be the upward neighbours of (x 1 , t) and (x 2 , t), respectively. For x 1 = x 2 , if we consider the layer above, since (x 1 , t) and (x 2 , t) have 3 upward neighbours each, the probability that (x 1 , t) or (x 2 , t) topples is equal to 1/3 times the number of topplings of their upward neighbours. Hence, Thus, from Eq. (9) we get for all Note that for the special casse of It is easy to verify that satisfies the linear equation (10), where G 2 [ X| Y ] denotes the expected number of topplings at X = (x, t), in an avalanche generated by a particle addition at Y = (y, t ). Clearly, The function To find a solution that is also consistent with the boundary conditions Eq. (12), we consider the following superposition of such solutions [2]: Here the summation over Y could extend over all sites. However, due to the light-cone structure of the propagator G 2 , only the sites which are in the future light-cone of O and in the past light-cones of both X 1 and X 2 contribute to the summation. Then, for X 1 = X 2 the equations determining the unknown coefficients f ( Y ) are the boundary conditions Eq. (12), which become where the summation over Y is over all sites that are in the future light-cone of O, and the backward light-cone of X. These equations are coupled linear equations which can be used to determine the unknown function f ( Y ) at sites Y in the past-light cone of X, and hence can be determined recursively starting from the driving site.
In fact, we do not need to know the f ( Y )'s for all Y ; the knowledge of their sum for each constant-time layer is sufficient. Let us denote this sum, for all times t ≥ 0, as Now, if we substitute our general solution from Eq. (15) into our expression for Φ 2 (t) as stated in Eq. (8) we get: Doing summations over x 1 , x 2 and t , using Eq. (7), and (17), we get Now let us define and note that K(t) vanishes with G 2 (x, t) when t < 0. Then, summing over different sites X in the layer t, in Eq. (16), and using Eq. (7) yields This equation differs from that in [2] by a factor, due to different choice of normalization of G 2 ( X) here. Also, note that the summation over t may be extended to +∞, as to be the generating functions of, K(z) and F (z), respectively. In terms of these generating functions, Eq. (21) can be expressed as Now let us remind ourselves that G 2 (x, t) is the expected number of topplings at site (x, t) given that there was a toppling at the origin O. A site can topple at most once, hence its expected number of topplings is exactly its probability to topple. Also, based on the Abelian property we know that the total number of topplings at another site (x , t ) is the sum of the topplings triggered by (x, t) toppling. Hence, we can write the expected number of topplings at site X = (x , t ) as for 0 ≤ t ≤ t . But then, the expected number of topplings at site (0, 2t) is just Substituting this in Eq. (20) yields Thus, the generating function becomes Now let us define H(z) = ∞ m=0 G 2 (0, m)z m and hence SinceK(z) is sum only over even values of t of G 2 (0, t), we haveK Substituting Eq. (5) with x = 0 into our definition of H(z) and evaluating the geometric sum yields Finally, evaluating this elementary integrals yields We can substitute this into, using Eq.(29), to get One can substituting this in Eq.(23), to getF (z) as an explicit function of z. Note that odd powers of √ z will cancel out in the expansion. From the fact that the dominant singularity ofK(z) for z near 1, is of the form (1 − z) −1/2 , we see that for z tending to 1 from below, the leading behavior ofF (z) also is (1 − z) −1/2 . Hence F (t) varies as t −1/2 for large t and then Φ 2 (t) varies as t 1/2 . This is the same behavior as found for the case K = 2 in [2]. The probability that Φ(t) is not zero decreases as t −1/2 for large t. But once it is non-zero, its typical value is of order t +1/2 , consistent with the mean value 1, see Eq. (7). Then the mean value Φ 2 (t) would be expected to grow as t 1/2 for large t. We now illustrate in Fig. 3 the importance of taking into account the corrections to scaling in estimating exponents from numerical data. First, we plot the exact values of Φ 2 (t) for 1 ≤ t ≤ 500. These values were determined using Eqs. (23) and (32) to expandF (z) as a Taylor series in z with Mathematica. A simple visual fit to a power-law gives Φ 2 (t) ≈ at α , with a ≈ 1.58 and α ≈ 0.52. Secondly, we plot the effective exponent α t . α t is defined in terms of the exact values of Φ 2 (t) at t and t + 1 by: We see that the effective exponent converges very slowly to the exact value 0.5. These two plots show that corrections to scaling are rather large in this problem. Conversely, knowing that Φ 2 (t) varies as t 1/2 implies that the probability that an avalanche has duration greater than t goes as t −1/2 . It follows that the avalanche duration exponent is 3/2. The avalanche dimension D in a directed model is identical to the duration exponent, that is D = 3/2 (Sec. 8.4.3 in [13]). All other exponents follow and we can verify that they are the same as the case of coordination number K = 2.
We note that a similar analysis has been reported for a directed sandpile model in [9]. In this paper, the authors determine the exact two-point correlation function G 2 (x, t) for a model with the following rules: the critical height is 4, and on toppling at (x, t), a particle is transferred to (x − 1, t + 1) and (x + 1, t + 1), and two particles to (x, t + 1). With these rules, the functions K(t) and F (t) were determined recursively numerically for small t (t < 500), and it was found that the corrections to scaling are large. However, no explicit analytical expressions for F (t) or K(t) were found.

V. NUMERICAL SIMULATIONS
The directed sandpile is comparatively easy to implement and study numerically even on very big lattices, because avalanches progress in one direction only, i.e., no backward avalanches [10] or multiple topplings occur [11]. The aim of the numerics is to provide reliable numerical estimates of (supposedly universal) critical exponents. In the following we will discuss pertinent issues in relation to the numerical simulations, the fitting models used and the critical exponent estimates found numerically.

A. Initialisation
The directed sandpile is deterministic up to the driving in the first row. As every state is recurrent, one might naively start from an empty lattice h(x, t) = 0 as this indeed belongs to the stationary state. Because every toppling moves K particles downstream and every particle added performs L moves before leaving the system, K particles additions cause L topplings (and thus KL moves) on average, i.e., the average avalanche size is where avalanche sizes s are measured as the number of topplings that occur in the system after a driving attempt (deposition of a particle in the top row). Note that this definition includes avalanche size s = 0. Starting from an empty initial configuration, the first moment of the cluster size s in a system with K = 4 and L = 512 (with M = 3072, see below) shows very good convergence within less than 2 · 10 7 avalanching attempts (i.e., particles deposited in the top row). However, the second moment of s still shows signs of drift after 5 · 10 10 avalanches. Higher moments show similar long transients, but are more noisy. This only underlines the fact that the steady state of model studied here has The data shown is for different K : 2 (red circles), 3 (green squares) and 4 (blue diamonds). For the same K, a higher curve corresponds to higher n. very slowly decaying temporal correlations, even though it has no spatial correlations.
To avoid these lengthy transients, we resorted to initialisation with random, independently, uniformly distributed h(x, t) < K. As the system is then no longer forced into an exceptional initial state, there is no noticeable drift in any of the moments measured. We have verified numerically for smaller lattices that random initial states generate the same moment estimates as starting from an empty system and taking estimates after those very long transients.

B. Parameters and Results
We used systems with sizes L = 2 r , with r taking integer values from 4 to 14, and K = 2, 3, 4 and widths M ≥ L(K −1) so that even the largest avalanche possible cannot topple all sites of a row. We used the Mersenne Twister pseudo-random number generator [12], which is (after initialising randomly the lattice) needed only to determine the site in the top row that receives a particle from the external drive. In order to determine scaling exponents, we measured moments s n of the avalanche size s for n = 2 to 5 (Fig. 4).
Statistical errors were determined by accumulating data over chunks of 10 6 avalanches and measuring the variance of those estimates. We produced usually at least several thousand chunks, except for the largest system sizes. If chunks are independent, then the variance does not vary noticeably when chunks are merged. Using that as an indicator, we merged chunks until they became independent. Based on the now independent measurements, statistical errors of the mean (across chunks) of a moment are given by the estimated standard deviation of the moments (across chunks), divided by the square root of the number of independent chunks.
Because high moments draw most of their weight from very large and thus very rare events, their relative error grows with their order. For example, for K = 3 and L = 16384, averaging over 3.8 · 10 8 avalanches, the fractional error in the eighth moment is approximately 0.014, while for the second moment it was approximately .00016. We report here our results only for moments up to order 5.
Standard finite-size scaling of the probability density of the avalanche sizes s implies that the moments scale to leading order in the system size L like [13] s n = a n L µn for L 1, where a n are metric factors, and µ n = D(1 − τ + n).
Here D is the fractal dimension of avalanche clusters, and τ ≥ 1 [14]. The critical exponents D and τ are expected to be universal, but they are not independent as s = L/K, Eq.(34) implies that D(2 − τ ) = 1.
For τ > 1 (which is expected in the present case) the scaling of the variance of the nth moment is dominated by that of the scaling of the 2nth moment and thus the relative error of independent samples scales like and thus grows, independent of the order, with the system size (height). Meaningful estimates of moments thus require larger and larger sample sizes for larger systems. However, simulation time per avalanche grows essentially like the average avalanche size, which is linear in the system size L. Worse, correlation times grow with the system size L, so that the increased demand on the sample size for larger system is met with highly increased costs for independent samples.
Using random initialisation, we were able to skip the transient as described above and thus could produce very large samples even for large system sizes. For the smaller system sizes, our sample sizes were comparatively large, typically about 10 10 avalanches and more. However, to fit well the accurate estimates for the moments for small L, we need to have a large number of corrections to scaling terms, which complicates the analysis. Therefore, we excluded system sizes smaller than L = 512 from further analysis.
First, the critical exponents D are extracted by applying moments analysis without including corrections to scaling. In the limit of large system sizes, the nth moment scales according to Eq. (35). Hence, plotting the measured moments s n against system size L yields estimates for µ n . According to Eq.(35) and the scaling relation τ = 2 − 1/D, we expect µ n = 1 + D(n − 1), so we can extract the critical exponent D by plotting µ n K D τ D (ctos) τ (ctos) 2 1.4999 (11)   I. Numerical estimates of the avalanche dimension D, and τ derived from it via the scaling relation τ = 2−1/D. The numerics for K = 2 were performed as a reference, as D = 3/2 and hence τ = 4/3 was already known analytically [2]. The first two columns are fits without including corrections to scaling. The last two columns list the exponents extracted when the fit includes corrections to scaling (ctos) Eq. (37) for systems with K = 3 and 4. The errorbars listed in the brackets correspond to three standard deviations of the last two significant digits.
vs. n. Then the avalanche size exponent is calculated form the scaling relation τ = 2 − 1/D. The critical exponents are listed in columns 1 and 2 in Tab. I. The critical exponents we find for K = 2 are consistent with the theoretical values D = 3/2 and τ = 4/3. However, the critical exponents for K = 3 and 4 show significant deviation from the theoretical values and we notice that the deviation increases with K.
However, we now include two corrections to scaling terms in our fitting model and fit each moment vs. system size L, using the Levenberg-Marquardt method [15]. Where convergence was a problem, we provided initial guesses by fitting first using fewer parameters (which typically results in substandard goodness of fit). The resulting estimates for exponents µ n were accepted for further analysis if the goodness of fit was at least 0.25 [15]. This was the case for n up to 5. One exponent, µ 1 , is known to be unity from Eq.(34), which means that µ n needs to be fitted linearly against 1 + D(n − 1) for n = 2, 3, 4, 5. Because the µ n are, for different n, all based on the same data (the set of avalanche sizes generated), they are bound to be correlated, but this is difficult to quantitfy reliably, except by using a rather brutal upper bound. Taking that course of action, we have effectively assumed that the estimates of µ n for n = 2, 3, 4, 5 may have been derived from distinct samples, by multiplying each error bar by a factor √ 4. The resulting estimates for the exponent D (and implicitly for τ ) are shown in Tab. I, columns 3 and 4, respectively.
These estimates fit acceptably well with the theoretical value of D = 3/2 and τ = 4/3, but only for very large sample sizes, system sizes and CPU time.

C. Summary and Conclusions
We have studied the directed Abelian sandpile model on a square lattice with K downward neighbors. When K = 2, avalanche clusters are compact without any holes. Using this property, the K = 2 case has previously been solved exactly [2]. When K > 2, avalanches clusters typically contain holes, that is, they are no longer compact. Hence, the previous derivation for the K = 2 case cannot be extended to the K > 2 case and it is an interesting question whether K > 2 belong to the same universality class as K = 2.
In this paper, we calculated exactly the exponents for the K = 3 case, where the problem is complicated by the fact that avalanche clusters are no longer compact. We find that the critical exponents are identical to the K = 2 case, that is, the avalanche dimension D = 3/2 and the avalanche size exponent τ = 4/3. In addition, we get exact expression for other observables, for example, generating functionF (z) of the mean-square flux. Our result shows that the deviations from K = 2 values observed in recent numerical studies are due to corrections to scaling.
We performed large scale numerical simulations of the generalized directed Abelian sandpile model for K = 2, 3 and 4. Although the empty state is a recurrent state, it is not a typical state. Initialising the system in the empty state results in extremely long transients before correlations caused by this exceptional state vanish. Initialising the system in a random recurrent state minimised the transient time before avalanche size moments estimates no longer drifted.
Using system sizes L = 512, 1024, . . . , 16384, moments analysis yields numerical estimates for the avalanche dimension D and hence τ = 2 − 1/D. For K = 2, the numerical estimates are consistent with the exact result. However, for K = 3 and 4, there are large corrections to scaling effects. If these are taken into account (see Eq.(37)), the resulting numerical estimates for the critical exponents are consistent with the exact findings, that is, the directed Abelian sandpile model belong to the same universality class for all K ≥ 2.