# Modified Renormalization Strategy for Sandpile Models

## Abstract

Following the Renormalization Group scheme recently developed by Pietronero et al, we introduce a simplifying strategy for the renormalization of the relaxation dynamics of sandpile models. In our scheme, five sub-cells at a generic scale $`b`$ form the renormalized cell at the next larger scale. Now the fixed point has a unique nonzero dynamical component that allows for a great simplification in the computation of the critical exponent $`z`$. The values obtained are in good agreement with both numerical and theoretical results previously reported.

preprint: HEP/123-qed  

The concept of Self-Organized Criticality (SOC) introduced by Bak, Tang and Wiesenfeld (BTW)  has attracted a wide interest to understand a class of dynamically driven systems which self-organize into a statistically stationary state characterized by the lack of any typical time or length scale. Numerical results for systems displaying SOC behavior have been extensively reported , but only a few theoretical approaches are known to be in agreement with numerical simulations in all dimensions. The major source of difficulties in the study of SOC systems lies in their inherent complexity that makes the models analytically tractable only in a few cases. The Abelian version of the BTW sandpile model, early addressed by Dhar , turned out to be one of these exceptions.

Recently, Pietronero, Vespignani, and Zapperi  developed a new type of real space renormalization group approach for dynamically driven systems, able to describe the self-organized critical state of sandpile models by defining a characterization of the phase space in which it is possible the renormalization of the dynamics under repeated change of scale. In addition, it is also possible to compute the critical exponents analytically . The method also reveals the nature of the SOC problems and provides a picture about the universality classes of different sandpile models. This scheme of renormalization has been recently improved by considering increasingly complex proliferation paths  and extended to forest-fire models .

In this report, we follow the renormalization procedure of references  but using a Greek cross-shaped cell instead of a square cell in the renormalization of the relaxation dynamics. The critical exponents that characterize the stationary state are then computed and they are found to be in good agreement with previous theoretical results and large scale numerical simulations both for the BTW and two state model of Manna. We will see that the use of this particular choice of cells simplifies the renormalization equations for the BTW model.

In what follows, we will focus on the sandpile critical height models in two dimensions. Sandpile models are cellular automaton defined on a lattice where to each site one assigns a variable (to which we will refer as energy). We let the system evolves by randomly adding units of energy on the system. When the energy of a site reaches a critical value, it relaxes releasing its entire energy to the neighboring sites. The affected sites may become unstable triggering new toppling events and so on until all sites are again stable. Three different classes of sites can be distinguished: (i) those sites for which the addition of a unit of energy does not induce relaxation (stable sites), (ii) those sites for which the addition of a unit of energy provokes they become unstable (critical sites), and (iii) unstable sites that will relax at the next time step. Open boundary conditions allow the energy to leave the system.

In this formalism, we will denote by $`\rho `$ the density of critical sites. These definitions can be extended to a generic scale $`b`$ by considering coarse grained variables. Thus, a cell at scale $`b`$ is considered critical if the addition of a unit of energy $`\delta E(b)`$ induces a relaxation of the size of the cell, that is, the subrelaxation processes span the cell and transfer energy to some neighbors. According to  the relaxation process can lead to four different possibilities at coarse grained levels: the energy can be distributed to one, two, three, or four neighbors with probabilities $`p_1`$, $`p_2`$, $`p_3`$, and $`p_4`$ respectively. Of course, it is also possible that in certain cases the unstable sites at coarse grained scale do not transfer energy to their nearest neighbors as well as to consider different proliferation problems. We, as in , will not consider these cases . Then, the probability distribution is defined by the vector

$$\stackrel{}{P}=(p_1,p_2,p_3,p_4)$$
(1)   
with the normalization condition $`_{i=1}^4p_i=1`$.

So, the properties of the system are fully characterized by the distribution $`(\rho ,\stackrel{}{P})`$ at each scale. The relation between $`\rho `$ and $`\stackrel{}{P}`$ can be derived by noting that in the stationary state the inflow of energy equals the flow of energy out of the system . This implies :

$$\rho ^{(k)}=\frac{1}{_iip_i^{(k)}},$$
(2)   
which allow us to evaluate the stationary distribution of critical sites at each scale $`k`$ of coarse graining.

Now, we define a renormalization transformation for the relaxation dynamics. We will use a cell-to-site transformation on a square lattice, in which each cell at scale $`b^{(k)}`$ is formed by five sub-cells at the finer scale $`b^{(k1)}`$ (see Fig. 1). We have chosen this type of cells for two reasons: one, because it implies the use of greater cells formed by five sub-cells at the finer scale, that is, when we scale up, five sub-cells form a new one at the larger scale; and second, one is intuitively tempted to follow the geometry of the relaxation that takes place in numerical simulations of sandpile models with energy transfer to N, E, S, W neighbors .

The length scaling factor is then $`\frac{b^{(k)}}{b^{(k1)}}=\sqrt{5}`$ (see Fig.1). Therefore, at a generic scale $`b^{(k)}`$, each cell is characterized by an index $`\alpha `$, ranging from one to five, indicating its number of critical sub-cells at the smaller scale $`b^{(k1)}`$. In order to ensure the connectivity properties of the avalanche in the renormalization procedure, only those configurations with three or more sub-cells at scale $`b^{(k1)}`$ can span the cell, transferring energy to $`i`$ neighboring cells. Thus, the starting relaxation processes $`p_i^{(k1)}`$ at scale $`b^{(k1)}`$ are renormalized in the correspondent process $`p_i^{(k)}`$ at scale $`b^{(k)}`$. Besides, it has been shown  that site correlations are averaged out in the stationary state. Therefore, taking into account this fact and the spanning rule we can write the weight of each configuration $`\alpha `$ in the stationary state as:

$`W_{(\alpha =3)}=2\rho ^3(1\rho )^2`$ (3)   
$`W_{(\alpha =4)}=4\rho ^4(1\rho )`$ (4)   
$`W_{(\alpha =5)}=\rho ^5`$ (5)   
Eq. (4) gives the probability that a cell at scale $`b^{(k)}`$ has the corresponding number of critical sub-cells at scale $`b^{(k1)}`$.

As an example of the general procedure, in Fig.2 we have drawn a series of relaxation processes $`p_1p_1p_2`$ at scale $`b^{(k1)}`$ that contributes to the renormalization of $`p_1^{(k)}`$ at the larger scale $`b^{(k)}`$, starting from a configuration of $`\alpha =3`$ critical sub-cells. The process consists of the following relaxation events that span the cell from left to right satisfying the spanning condition. First, the unstable sub-cell on the left relaxes toward the other critical sub-cell (the center one, Fig.2b) which occurs with probability $`(1/4)p_1^{(k1)}`$, where the index $`(k1)`$ denotes that the relaxation takes place at scale $`b^{(k1)}`$. Second, we consider the process in which the new unstable sub-cell also relaxes toward the sub-cell on the right through another $`p_1`$ process (Fig.2c) which again happens with a probability $`(1/4)p_1^{(k1)}`$. Finally, the sub-cell on the right has become unstable and transfers with probability $`(2/3)p_2^{(k1)}`$ two units of energy one inside and one outside the original cell of size $`b^{(k)}`$ (Fig.2d). The series of processes described contributes to the renormalization of $`p_1^{(k)}`$. Nevertheless, it is necessary to note that the relaxations displayed in Fig 2a-2d are not all the processes that contribute to the renormalization of $`p_1^{(k)}`$ through a $`p_1p_1p_2`$ series. Fig.2e shows a $`p_2`$ relaxation event that, although involves two neighboring sites outside the original cell of size $`b^{(k)}`$, also contributes to the renormalization of $`p_1^{(k)}`$ with probability $`(1/6)p_2^{(k1)}`$. This is a new characteristic inherent to the cell-to-site transformation chosen. Now, if we take into account all the processes that lead to $`p_1^{(k)}`$, for $`\alpha =3`$, one gets

$`p_1^{(k)}={\displaystyle \frac{1}{3}}\left\{\left({\displaystyle \frac{1}{6}}p_2^{(k1)}+{\displaystyle \frac{1}{2}}p_3^{(k1)}+p_4^{(k1)}\right)\left({\displaystyle \frac{1}{4}}p_1^{(k1)}\right)\left({\displaystyle \frac{3}{2}}p_1^{(k1)}+{\displaystyle \frac{4}{3}}p_2^{(k1)}+{\displaystyle \frac{1}{2}}p_3^{(k1)}\right)\right\}`$ (6)   
$`+{\displaystyle \frac{2}{3}}\left\{\left({\displaystyle \frac{1}{4}}p_1^{(k1)}+{\displaystyle \frac{1}{2}}p_2^{(k1)}+{\displaystyle \frac{3}{4}}p_3^{(k1)}+p_4^{(k1)}\right)\left({\displaystyle \frac{1}{4}}p_1^{(k1)}\right)\left({\displaystyle \frac{3}{4}}p_1^{(k1)}+{\displaystyle \frac{7}{6}}p_2^{(k1)}+{\displaystyle \frac{1}{2}}p_3^{(k1)}\right)\right\}`$ (7)   
where in Eq. (7) the factors $`\frac{1}{3}`$ and $`\frac{2}{3}`$ refer to the multiplicities of the configurations (see Fig.3).

In a similar way (though much more complicated), one obtains expressions for $`p_i^{(k)}`$, $`(i=2,3,4)`$, for $`\alpha =3`$ and imposes the normalization condition $`_{i=1}^4p_i^k=1`$. The procedure is repeated taking into account the configurations with $`\alpha =4`$ and $`\alpha =5`$ critical sites and the renormalized probabilities at level $`k`$ are finally derived by averaging over the configurations of different $`\alpha `$-values including their statistical weights $`W_\alpha (\rho ^{(k1)})`$. Therefore, the probabilities $`p_i^{(k)}`$ at length scale $`b^{(k)}`$ will be given by

$$p_i^{(k)}=\underset{\alpha =3}{\overset{5}{}}W_\alpha (\rho ^{(k1)})p_i^{(k1)}(\alpha )$$
(8)   
with $`W_\alpha (\rho ^{(k1)})`$ and $`\rho ^{(k1)}`$ given by Eq. (4) and Eq. (2), respectively. As the computation of the probabilities $`p_i^{(k)}`$ in Eq. (8) is rather lengthy and cumbersome, we have developed a C-code to compute all the polynomial term coefficients that contribute to the renormalization transformation.

Now, we proceed to explore the scale-invariant behavior of the model by finding the fixed-point solution $`p_i^{(k1)}=p_i^{(k)}`$. In order to do this, we start from the shortest length scale characterized by $`(\rho ^{(0)},\stackrel{}{p^{(0)}})`$ and study how it evolves under repeated iteration of the transformation (8). For the two state model of Manna  the parameters $`(\rho ^{(0)},\stackrel{}{p^{(0)}})`$ are $`(\rho ^{(0)},0,1,0,0)`$ whereas for the BTW sandpile we have $`(\rho ^{(0)},0,0,0,1)`$. Here, the initial value of the density of critical sites $`\rho ^{(0)}`$ is irrelevant for the dynamics since the system evolves to a fixed point whatever be the distribution of critical sites at the small scale dynamics.

As in Refs , both models have the same fixed point, but here there is an important difference in relation to the value of the fixed point. We obtain for the fixed point the value $`(\rho ^{},\stackrel{}{p}^{})=(\frac{1}{4},0,0,0,1)`$, that is, in the BTW model one starts from the fixed point! This is not indeed the case for the two-state model of Manna for which we need to iterate Eq. (8) more than twenty times to reach the same fixed point. We believe that this is a consequence of our renormalization strategy for the relaxation dynamics and constitutes a great simplification in the calculation of the dynamical exponent $`z`$. In fact, we were expecting the existence of a critical fixed-point value different from that reported in references  although the critical exponents should be very close since they are determined by the properties of the system at large scales.

The exponent $`\tau `$ that characterizes the power-law avalanche size distribution can be obtained following the procedure of . Consider the probability $`K_{b^{(k1)},b^{(k)}}`$ that the relaxation processes that are active at scale $`b^{(k1)}`$ do not extend beyond the larger scale $`b^{(k)}`$. This is expressed as 

$$K=\frac{\underset{b^{(k1)}}{\overset{b^{(k)}}{}}P(r)𝑑r}{\underset{b^{(k1)}}{\overset{\mathrm{}}{}}P(r)𝑑r}=1\left(\frac{b^{(k)}}{b^{(k1)}}\right)^{2(1\tau )}=1\left(\sqrt{5}\right)^{2(1\tau )}.$$
(9)   
Eq. (9) also satisfies

$$K=p_1^{}(1\rho ^{})+p_2^{}(1\rho ^{})^2+p_3^{}(1\rho ^{})^3+p_4^{}(1\rho ^{})^4$$
(10)   
Then, the exponent $`\tau `$ is given by

$$\tau =1\frac{1}{2}\frac{\mathrm{log}(1K)}{\mathrm{log}(\sqrt{5})}=1.235.$$
(11)   
This value of $`\tau `$ is in very good agreement with the value obtained in  and with large-scale numerical simulations which give $`\tau =1.27`$ for the two-state model of Manna, and $`\tau =1.29`$ for the BTW sandpile model .

A second independent critical exponent can also be computed. This is the so-called dynamical exponent $`z`$ that relates the spatial scale $`r`$ to time scale $`t`$ through the power law $`tr^z`$. As pointed out in , the calculation of $`z`$ could be an enormous and laborious task because the knowledge of the fixed point value is not sufficient and we have to know the complete form of the renormalized dynamics. Nevertheless, as we said before, the use of our larger cells in the renormalization transformation leads to a fixed point with a unique nonzero component in the vector $`\stackrel{}{p}^{}`$. This constitutes a great simplification in the derivation of the complete structure of the renormalized dynamics. In what follows, we will derive at a glance the dynamical critical exponent for the BTW sandpile model. In order to obtain the dynamical exponent we have to calculate the average number $`<t>`$ of noncontemporary processes at scale $`b^{(k1)}`$ needed to have a relaxation process at the larger scale $`b^{(k)}`$, which is related with $`z`$ through

$$z=\frac{\mathrm{ln}<t>}{\mathrm{ln}\left(\frac{b^{(k)}}{b^{(k1)}}\right)}=\frac{\mathrm{ln}<t>}{\mathrm{ln}\left(\sqrt{5}\right)}.$$
(12)   
In Fig.3 we have depicted the possible starting configurations for the different values of $`\alpha `$. It is also shown the time steps needed to have a relaxation process at the larger scale. Such a simplification in the calculus is possible because we have to consider only the relaxations that contribute to the renormalization of $`p_4`$ at larger scale. As can be seen, we need two time steps for the symmetric configurations (those in which the initial unstable site is located at the center of the cell) and three for the non-symmetric configurations (those in which the initial unstable site is located in one of the critical boundary sites of the cell). Therefore,

$$<t>=\frac{1}{\underset{\alpha }{}W_\alpha (\rho )}\underset{\alpha }{}t^{}(\alpha )W_\alpha (\rho )$$
(13)   
where $`t^{}(\alpha )`$ is the weighted average of time steps taking into account the different additional statistical weights due to multiplicities $`\omega `$ in each configuration $`\alpha `$ (see Fig.3). Now, evaluating Eq. (13) at the fixed point we obtain,

$$z=1.236$$
(14)   
The value (14) is in remarkably good agreement with the numerical result $`z=1.21`$  and with the exact value $`z=5/4`$ . The other critical exponents can be derived from scaling relations . Table I summarizes the values of the critical exponents obtained for the BTW sandpile model and those reported by previous renormalization scheme and numerical simulations.

In this report, we have introduced an alternative renormalization strategy that simplifies the analytical derivation of the critical exponents that characterize the dynamics of sandpile models. By using larger cells, formed by five sub-cells of the finer scale, we obtain a fixed point with a unique nonzero dynamical component which allow us to derive the whole form of the renormalized dynamics in a more direct and simple way. The values of the exponents here obtained are in good agreement with those previously reported. Besides, as in similar analytical predictions, the two-state model of Manna and the BTW sandpile model belong to the same universality class . The results confirm the robustness of the renormalization group approach.

It is a pleasure to thank A. Vespignani for stimulating discussions and G. Caldarelli for useful correspondence. Y.M would like to thank the AECI for financial support.