# Fluctuations and correlations in sandpile models

## Abstract

We perform numerical simulations of the sandpile model for non-vanishing driving fields $`h`$ and dissipation rates $`ϵ`$. Unlike simulations performed in the slow driving limit, the unique time scale present in our system allows us to measure unambiguously response and correlation functions. We discuss the dynamic scaling of the model and show that fluctuation-dissipation relations are not obeyed in this system.

PACS numbers: 05.65.+b, 05.70.Ln

The sandpile automaton  is one of the simplest model of avalanche transport, a phenomenon of growing experimental and theoretical interest. In the model introduced by Bak, Tang and Wiesenfeld (BTW) , grains of “energy” are injected into the system. Open boundary conditions  or bulk dissipation insure a balance between input and output flow and allow for a non-equilibrium stationary state. In the limit of slow external driving and small dissipation, which corresponds to an infinite time scale separation between driving and response, the model displays an highly fluctuating avalanche behavior, indicative of a critical point. Despite the impressive theoretical effort devoted to understanding the critical behavior of the model , several important issues still remain to be addressed.

Numerical simulations are usually performed under slow driving and boundary dissipation, since the limit of infinite time scale separation is easily implemented in the computer and provides a simple way to access the avalanche critical behavior . However, due to the presence of two infinitely separated time scales, an unambiguous definition of dynamic response and correlation functions is not possible . This hinders a clear characterization of the non-equilibrium stationary state in terms of static and dynamic response and correlation functions. Evaluation of these quantities helps to elucidate the nature of the critical point and provides a test of fluctuation-dissipation relations, at least in some weaker sense. Recently, it has been proposed to interpret the behavior of sandpile models in analogy with other non-equilibrium critical phenomena, such as absorbing phase transitions , driven interfaces in random media  and branching processes . These theoretical studies suggest new ways to perform numerical simulations in which a unique time scale is considered .

In this letter, we present numerical simulations of the sandpile model for different driving rates $`h`$ and study how the system approaches the critical point when $`h0`$. In this way, we are able to measure quantities that are not accessible in the time scale separation regime. The local density of active sites, that can be identified as the order parameter of the model , is homogeneous only in the case of bulk dissipation. For boundary dissipation, it displays a marked curvature, that was anticipated in Refs. and could explain several scaling anomalies found in the BTW model. The energy landscape is instead homogeneous in both cases and its statistical properties do not depend on the dissipation rate $`ϵ`$ in the limit $`h0`$.

We measure correlation and response functions in time and space domains and observe the scaling of the related characteristic lengths and times. We find two different characteristic times, implying that fluctuation-dissipation relations are not obeyed. We observe, however, a well defined scaling behavior and the value of the critical exponents are in agreement with recent large scale numerical simulations of slowly driven sandpiles . Finally, the present numerical analysis opens the road to future studies to resolve some longstanding problems such as the precise identification of universality classes for these models .

In sandpile models , each site $`i`$ of a $`d`$dimensional lattice bears an integer variable $`z_i0`$, which we call energy. At each time step an energy grain is added on a randomly chosen site. When a site reaches or exceeds a threshold $`z_c`$ it topples: $`z_iz_iz_c`$, and $`z_jz_j+1`$ at each of the $`g`$ nearest neighbors (nn) of $`i`$. Each toppling can trigger nn to topple and so on, generating an avalanche. The original BTW model is conservative and energy is dissipated only at the boundary, i.e. energy grains from toppling boundary sites flow out of the system. Infinitely slow driving is implicitly built into the model: during the avalanche the energy input stops, until the system is again quiescent (no active sites are present), so that we can identify two distinct time scales $`T_d`$ and $`T_a`$, for driving and activity, respectively. A single driving time step can in principle be followed by an infinite number of avalanche time steps and $`T_a/T_d0`$. For this reason, there are two possible definitions for the correlation function, depending on the choice of the scale used to measure time (slow or fast) .

Here we simulate the BTW sandpile model for a non-vanishing driving field: each site has a probability $`h`$ per unit time to receive an energy grain, also if active sites are present in the system. This defines a unique time step for both driving and activity updating. The parameter $`h`$ sets the driving rate, and in the limit $`h0^+`$ we recover the slow driving limit; i.e. during the evolution of an avalanche the system does not receive energy. We consider two possible mechanisms for energy dissipation: (i) usual boundary dissipation and (ii) bulk dissipation, simulated introducing the probability $`\alpha `$ that a toppling site looses its energy without transferring it to the neighbors, which corresponds to an effective average dissipation $`ϵ=\alpha z_c`$. In case (ii), we impose periodic boundary conditions. We use two dimensional lattices with linear sizes ranging from $`L=64`$ to $`L=300`$, and parameters in the range $`10^6<h<10^3`$ and $`10^3<ϵ<10^2`$.

The order parameter in sandpile models is the density $`\rho _a`$  of active sites, whose energy is larger than $`z_c`$ . The dependence of the order parameter on the control parameters $`h`$ and $`ϵ`$ is readily obtained by means of conservation arguments : since energy is conserved in the stationary state, the incoming energy flux $`J_{in}=hL^d`$ must be equal to the dissipated energy $`J_{out}=ϵ\rho _aL^d`$. By equating the two fluxes we obtain $`\rho _a=h/ϵ`$. In systems with boundary dissipation, the effective dissipation rate scales with the system size as $`ϵL^\mu `$, with $`\mu =2`$ , yielding $`\rho _ahL^2`$. It has to be noticed that the model is critical just in the double limit $`ϵ0`$ and $`h/ϵ0`$. The onset of a nonvanishing driving thus destroys criticality in that it enforces a nonzero dissipation. For $`hϵ`$ the cutoff length scaling is dominated only by dissipation, while for greater driving fields more complicate scaling behaviors occur.

We study the behavior of the density of active sites in the system and measure the stationary average density of local energies; i.e. the density $`\rho _i`$ of sites with $`i`$ energy grains. In Fig. 1 we report the behavior of the densities as a function of $`h/ϵ`$. For small values of $`h/ϵ`$ we find $`\rho _i=\rho _i^0+c_ih/ϵ+𝒪((h/ϵ)^2)`$, where $`\rho _i^0`$ are the values extrapolated from the limit $`h0^+`$ and are given by $`\rho _0^0=0.075(1)`$, $`\rho _1^0=0.176(1)`$, $`\rho _2^0=0.307(1)`$ and $`\rho _3^0=0.442(1)`$. These values are in excellent agreement with the exact values obtained for the slowly driven sandpile (with boundary dissipation)  and are independent of the dissipation rate. This implies that the energy substrate over which avalanches propagate is the same in the case of bulk and boundary dissipation. For $`i>3`$ we obtain $`\rho _i^0=0`$ and for small $`h`$ we observe $`\rho _a\rho _4`$, while for larger $`h`$ higher energy levels become populated and $`\rho _a`$ has non negligible contributions coming from $`\rho _i`$ with $`i>4`$. Finally we confirm that $`\rho _a=h/ϵ`$ for the whole range of parameters. In the case of boundary dissipation we recover $`\rho _ahL^2`$.

To elucidate the differences between bulk and boundary dissipation, we measure the local density of active sites $`\rho _a(r)`$. In the case of bulk dissipation the density profile is flat $`\rho _a(r)=\rho _a`$, while in the case of boundary dissipation we obtain a surface that can be well approximated by a paraboloid (see Fig. 2). This is due to the highly inhomogeneous dissipation which imposes a zero density of active sites on the lattice boundary and corresponds to an elastic interface pinned at the boundaries as discussed in Ref.s . This effect can explain the anomalies encountered in the numerical evaluation of avalanche exponents  and the persistent deviations from simple scaling observed in BTW model .

In order to obtain a quantitative description of the stationary state, we study the effect on the stationary density of a small perturbation in the driving field

$$\mathrm{\Delta }\rho _a(r,t)=\chi _{h,ϵ}(rr^{},tt^{})\mathrm{\Delta }h(r^{},t^{})𝑑r^{}𝑑t^{}$$
(1)   
where $`\chi _{h,ϵ}(r,t)`$ is the local response function. In the limit $`h0^+`$ the integrated susceptibility $`\chi 𝑑td^dr\chi _{h,ϵ}(r,t)`$ scales as the average avalanche size $`\chi s`$ and the time integrated response function scales as 

$$\overline{\chi }_{h0,ϵ}(r)\chi _{h0,ϵ}(r,t)𝑑t\frac{1}{r^{d2}}e^{r/\xi }$$
(2)   
where $`\xi `$ is the characteristic length. Since $`\chi =\rho _a/h`$ and $`\rho _a=h/ϵ`$, the response function diverges in the limit of vanishing driving and dissipation as $`\chi =1/ϵ`$. By noting that $`\chi \xi ^2`$, we obtain that $`\xi ϵ^\nu `$ with $`\nu =1/2`$ . These results can be obtained in mean-field (MF) theory but hold in all dimensions due to conservation .

To measure the response function, we drive the system in the stationary state with a given $`h`$ and we then add $`n`$ energy grains (i.e. $`\mathrm{\Delta }h=n/L^2`$) on a given lattice site. The time integrated response function is equivalent to the average difference of activity $`\overline{\chi }_{h,ϵ}(r)=\mathrm{\Delta }\rho _a(r)=\rho _a(r)_{h+\mathrm{\Delta }h}\rho _a(r)_h`$, where $`r`$ denotes the distance from the perturbed site. We observe that this function decays exponentially as predicted by Eq. 2 and measure the correlation length $`\xi `$ (see Fig. 3). In the case of bulk dissipation, for small driving fields the $`\xi `$ depends only on the dissipation rate and scales as $`\xi ϵ^\nu `$, with $`\nu =0.50\pm 0.01`$ (see Fig. 4). In the case of boundary dissipation, to evaluate $`\overline{\chi }_{h,L}(r)`$ we have to consider explicitly the spatial inhomogeneity of the stationary density: $`\rho _a(r)h/ϵ`$. We observe also in this case that the integrated response function decays exponentially and defines a correlation length increasing linearly with the lattice size; i.e $`\xi L`$. This result does not agree with the anomalous scaling found in a continuous energy sandpile. We perform analogous simulations in $`d=3`$ and find that Eq. 2 is still verified .

Furthermore, we study the response function in the time domain defined as $`\stackrel{~}{\chi }_{h,ϵ}(t)𝑑r\chi _{h,ϵ}(r,t)`$ after a small variation $`\mathrm{\Delta }h`$ of the driving field. Also in this case we obtain a clear exponential behavior defining the characteristic time scale $`\tau `$. For small driving field, $`\tau `$ scales as a function of the dissipation rate as $`\tau ϵ^\mathrm{\Delta }`$, with $`\mathrm{\Delta }=0.75\pm 0.05`$ (see Fig. 4). We then evaluate the dynamical exponent $`z=\mathrm{\Delta }/\nu =1.5\pm 0.1`$ relating time and spatial characteristic length: $`\tau \xi ^z`$. In the limit $`h0^+`$, we expect that the critical exponents $`\nu `$ and $`z`$ express the divergence of avalanche characteristic size and time, respectively. The numerical results confirm this observation . For increasing values of $`h`$, the driving field enters the scaling form and the results will be reported elsewhere.

We now turn to the analysis of the correlation function defined as $`C(r,t)=\rho _a(r,t)\rho _a(0,0)\rho _a^2`$. In previous simulations, performed in the slow driving limit, correlation functions were usually measured with respect to the slow time scale  and the fast time scale was explored studying the avalanche propagation. The introduction of non vanishing driving and dissipation allows us to bridge the gap between the two regimes. We study the correlation function in time and space domains and find an exponential decay at long times and distances , defining the correlation lengths $`\xi _c`$ and $`\tau _c`$ for space and time, respectively. The scaling of these correlation lengths is in agreement with the one obtained analyzing the response functions (i.e. $`\xi _cϵ^\nu `$, with $`\nu 0.5`$ and $`\tau _cϵ^\mathrm{\Delta }`$, with $`\mathrm{\Delta }0.75`$) and confirms the existence of a unique critical behavior in time and space (see Fig. 4).

In order to clarify the interplay between slow and fast dynamical modes, we analyze fluctuation-dissipation relations. In equilibrium phenomena, the fluctuation-dissipation theorem ensures that the response of the system to a small perturbation is related to the correlation function. In particular, the response function is given by

$$\chi (t)=\frac{1}{T}\frac{dC(t)}{dt},$$
(3)   
where $`T`$ is the temperature. Eq. 3 is strictly verified only in equilibrium systems, but it has been recently generalized to some classes of non equilibrium systems, namely systems displaying “aging” . In those examples the fluctuation-dissipation relation provides an information on an effective non equilibrium temperature that rules the dynamical evolution of the system.

We test Eq. 3 and we find that the usual linear behavior does not hold. On the contrary, we show that the parametric plot of $`\chi (t)`$ versus $`C(t)`$ defines a power law behavior, as shown in the double logarithmic plot of Fig. 5. This is striking evidence that the fluctuation-dissipation relation does not hold in these systems. Since we are in presence of two exponential functions, the linear behavior on the logarithmic scale is the signature of two different values for characteristic times $`\tau _c`$ and $`\tau `$ for the correlation and response function respectively. The slope indicates the ratio among the two time scales is given by $`\tau _c/\tau 0.4`$ and does not depend on driving and dissipation rates. This observation reflects the fact that the correlation and the response times scale with the same exponents with respect to dissipation and define unambiguously the critical behavior of the model. In particular, it implies that the dynamical exponent $`z1.5`$ is unique and can be estimated either by measuring avalanche distributions or the correlation functions. Previous simulations  revealed two different dynamical exponents in the fast and in the slow time scale. These differences are probably due to the ambiguous definition of time in the infinite time scale separation limit.

Finally, we note that it is not possible to define an effective temperature for the dynamics of sandpile models. It is interesting to compare this observation with a recent work  showing that the stationary state of non-equilibrium threshold models, similar to the one studied here, can be described by Boltzmann statistics in the mean-field limit. The validity of claim of Ref.  has been debated in the literature . We measure fluctuation-dissipation relations in a random neighbor sandpile model, which is described by mean-field theory, and find that fluctuation-dissipation relations are not satisfied , in disagreement with the conclusions of Ref. .

We thank A.Chessa, D.Dhar, R. Dickman, S. Franz, K.B. Lauritsen, E.Marinari, M. A. Muñoz, R. Pastor-Satorras and A. Stella for comments and discussions. A.V. and S.Z. acknowledge partial support from the European Network Contract ERBFMRXCT980183.