Time-Dependent Properties of Sandpiles

Bak, Tang, and Wiesenfeld (BTW) proposed the theory of self-organized criticality (SOC), and sandpile models, to connect “1/f” noise, observed in systems in a diverse natural setting, to the fractal spatial structure. We review some of the existing works on the problem of characterizing time-dependent properties of sandpiles and try to explore if the BTW's original ambition has really been fulfilled. We discuss the exact hydrodynamic structure in a class of conserved stochastic sandpiles, undergoing a non-equilibrium absorbing phase transition. We illustrate how the hydrodynamic framework can be used to capture long-ranged spatio-temporal correlations in terms of large-scale transport and relaxation properties of the systems. We particularly emphasize certain interesting aspects of sandpiles—the transport instabilities, which emerge through the threshold-activated nature of the dynamics in the systems. We also point out some open issues at the end.


INTRODUCTION
More than three decades ago, Bak, Tang, and Wiesenfeld (BTW) proposed the theory of "selforganized criticality" (SOC) [1], and sandpile models [2,3], as an explanation of the physical origin of the spatio-temporal scale invariance in natural systems found around us [4][5][6]. Indeed, in nature scale invariance is rather a rule, than an exception [7,8]. One abundantly finds fractal spatial structures, such as mountain ranges, coastlines and river basins [9,10], etc., which can be characterized by power-law correlations. For example, consider the height-height correlation function [h(x) − h(x ′ )] 2 ∼ |x − x ′ | 2χ in a mountain range, where h(x) and h(x ′ ) are heights at position x and x ′ . Clearly, the correlation function varies with the distance |x − x ′ | as a power law, characterized by an exponent 0 < χ < 1. In other words, upon varying the length scale, mountain ranges would appear a "self-similar" object: if one zooms in or out, a picture of a mountain would look almost the same. Likewise, the bird's eye view of the structures at the delta area of a river, such as the Ganges can look quite similar to the eyes if one changes the length scales (simply by broadening the horizon); similarly, the satellite picture of a country's coastline could well appear self-similar at different length scales. Indeed, such self-similarities indicate that these natural objects generate long-ranged power-law correlations in space.
Similar nontrivial structures can be observed also in the time (or, equivalently, frequency) domain, where many natural phenomena exhibit scale-invariant behaviors, implying the existence of long-ranged temporal correlations in the systems. For example, consider voltage fluctuations across a resistor in a current-carrying conductor, operating in a non-equilibrium steady-state condition. Due to the inherent noise in the system, the voltage V(t) across the resistor randomly fluctuates as a function of time t, however has a steady value V(t) = V 0 . Now let us consider the voltage power spectrum, i.e., the Fourier transform S(f ) ∼ e −i2π ft C(t)dt of the voltage autocorrelation function C(t) = V(t)V(0) − V 2 0 . Indeed, in a broad range of frequencies (10 −7 Hz < f < 10 4 Hz) and independent of any specific nature of the systems, the power spectrum has a power-law form as S(f ) ∼ 1/f α with the exponent generally lying in the range 1 < ∼ α < ∼ 1.4; for a review of the 1/f noise in solids (see [11]). This low-frequency noise is called in the literature the "flicker" or, simply, the "1/f " noise, where the exponent can in principle be in the range 0 < α < 2 (the case with α = 0 corresponds to the white or delta-correlated noise; on the other hand, α < 0 implies anti-correlation in a time-signal). Notably, the 1/f noise is different from the Johnson-Nyquist noise [12], which does not depend on frequency (α = 0 in that case) and occurs in equilibrium systems; it is different from the shot noise, which is also independent of frequency and happens near equilibrium. While it is less erratic than the white noise, the 1/f noise is somewhat more erratic than the 1/f 2 noise. Quite remarkably, the 1/f noise has been observed in a diverse range of apparently unrelated non-equilibrium phenomena, such as in traffic movements [13], flow of sand in hour glasses [14], and solar flare [15], etc. Not surprisingly, the origin of the subtle longtime correlations in 1/f noise has lacked a general theoretical understanding so far.
In their original paper [1], BTW attempted to provide a universal mechanism of the 1/f noise through the interplay between scale-invariant spatial and temporal structures, which can develop in slowly driven spatially extended dynamical systems observed in nature. There were two key ingredients in their theory: threshold-activated dynamics and slow driving. Bak, Tang, and Wiesenfeld argued that these systems, being highly non-linear, are governed by dynamics where there are dynamical activities only when the value of a certain local quantity crosses a threshold value; the system dissipates energy to maintain itself in a non-equilibrium steady state. According to BTW, under slow driving, such systems would "self-organize" to a "minimally stable" state, which is highly sensitive to perturbations and can create fast relaxation activities, called avalanches. The activities can occur at all length scales, implying a critical state for the system. Here the "slow" driving essentially refers to the separation of time scales in the systems: on the time scale of the driving, the systems relax instantaneously through avalanches. BTW dubbed this particular mechanism of achieving a critical state as the "self-organized criticality" (SOC). That is, apparently it does not require fine-tuning of any external parameter(s) to maintain the criticality; this is unlike any equilibrium critical phenomena, which are achieved only by tuning, say, temperature, and chemical potential (or, magnetic field in a magnetic system).
To demonstrate the mechanism of SOC, BTW proposed a model, called the BTW sandpile [1,16], which is nothing but simply a metaphor for a real sandpile, made on a finite base and driven to the edge of its stability through slow addition of grains and dissipation (loss of grains) at the boundary. On a two dimensional square lattice, the model is defined as follows. At any site (x, y), a height variable h(x, y, t) at time t is defined as the number of grains at the site. If h(x, y, t) ≥ h c = 2d (dimension d = 2 here), 2d grains topple and each of the neighboring sites get one grain each, i.e., h(x, y ± 1, t + 1) = h(x, y ± 1, t) + 1. (1) All sites are updated in parallel and the dissipation is incorporated through the boundary condition by keeping h(x, y, t) = 0 at the boundary. Note that the number of grains are conserved in the bulk; grains are dissipated only at the boundary.
As the toppling condition depends on the height variable, in this article we shall refer to this class of models as critical height-type one. The model can be straightforwardly generalized to higher (and lower) dimension. Later several stochastic variants of the BTW model, such as the two celebrated models-the Manna sandpile [17] and the Oslo ricepile [18], with stochastic toppling rules, were introduced. It has been debated whether the SOC systems can be thought of as ones spontaneously evolving, or "self-organizing." toward criticality, especially when the driving rate itself is "tuned" to zero [19]. (Nevertheless, one could simply take a view that such slowly driven systems exist in nature and are worth exploring for their interesting properties.) To elucidate this issue of "selforganization" vs. "fine tuning", the fixed-energy variants of the BTW sandpiles were introduced [20][21][22]. In the fixed-energy sandpiles, there is no dissipation of grains and the total mass in the systems remains conserved. Upon tuning the density below a critical density, the system is observed to undergo an absorbing phase transition (APT) [23], where all dynamical activities in the system ceases to exist. Indeed, in the stochastic fixed-energy sandpiles, the critical density is found to be the same as the steady-state density in the corresponding slowly driven dissipative sandpiles [24]; however, the situation can be complicated in the case of deterministic sandpiles, such as the BTW model, due to the ergodicity breaking [25,26].
Despite the fact that the original motivation of BTW was to connect the long-ranged spatial and temporal correlations observed in natural phenomena, much attention has been drawn in the past toward characterizing the spatial structures of sandpiles, mainly through the avalanches statistics [27]. Indeed, there is a vast body of experimental [28,29], simulation [17,21,22,[30][31][32][33][34][35][36][37][38][39][40][41][42] and exact [16,[43][44][45][46][47][48][49][50] results available in the literature; for reviews, see [19,51]. The time-dependent properties, on the other hand, have not been explored as much. In this article, we review some of the past works, which attempted, and achieved to a certain extent, dynamical characterization of sandpiles, through studies of simple models, which are amenable to analytic calculations and easy to simulate on computers. However, the task of large-scale hydrodynamic characterization of relaxation and transport properties of sandpiles is still incomplete. This review would perhaps give some hints of the possible directions in which one could proceed.
Organization of the article is as follows. In section 2, we discuss the nature of the long-range temporal correlations present in sandpiles, through studies of power spectrum of various time-dependent quantities, which have been measured in the experiments and simulations in the past. In this context, we also discuss the theories, which have been developed to capture these long-time correlations in the system. In section 3, we review some of the other interesting aspects of timedependent correlations in the context of particle transport, i.e., that of tagged-particle diffusion in sandpiles. We then address the large-scale relaxation and the transport instabilities in sandpiles through an exact hydrodynamical framework. We end the article with concluding remarks and some open issues.

Experiments
After the concept of SOC was introduced, it was quite natural to ask whether a pile of real sand, or any other granular material for that matter, exhibits SOC or not. This question was addressed by Jaeger et al. in an experiment [28] on a pile of granular materials (glass beads and aluminum-oxides grains) driven in a slowly rotating drum, which tilts the surface of the pile and thus generates avalanches of grains down the slope of the pile. The outflux of the grains (i.e., particle-current at the boundary) as a function of time, and the power-spectrum of the corresponding time-signal, were calculated by measuring the change of capacitance (which is related to the particle flow-rate) in a parallel-plate capacitor through which the particles passed. The power-spectrum was observed to have a broadened peak with a subsequent 1/f 3 decay for large frequencies; the peak in the power-spectrum corresponds to the average interval between two successive avalanches and the large-frequency power-law decay essentially arose from the short-time correlations in the particle currents. In other words, the non-trivial 1/f α spectrum with 1 < ∼ α < 2 was not observed in the experiment. For a real granular pile, the 1/f noise is perhaps not expected to occur, due to the dominant effects of inertia of grains, which tend to roll down the slope under gravity. In the presence of dissipation, the inertia essentially introduces a crossover length scale, beyond which the dissipative regime takes over. This cross over length scale could however be large in a typical experimental set up. To reduce the inertial effects, the Oslo group [29] studied avalanches in a pile of long-grained rice in a specially designed experiment. The pile was made between two vertical plates kept on a plane base, where the separation between the plates is small compared to the length of the base. One end of the pile is closed and the other end is open; the grains are added slowly through the closed end and the grains go out of the system through the open end. Large anisotropy of the grains stops the grains from rolling too far down the surface of the pile. In this particular set up, SOC was indeed observed in the avalanche statistics, where avalanche sizes were observed to follow a power-law distribution, implying events at all length scales and hence criticality. Here avalanche size is defined as the spatial extent of activities generated in the system when the system is perturbed by addition of a single grain (in experiment, the avalanche size can be measured, e.g., through the number of displaced grains in response to a perturbation localized in space and time). The ricepile experiment shows that SOC may not be that generic a feature of systems having a threshold activated dynamics, rather depends on various other details (e.g., shapes and sizes of grains and what materials they are made of, etc.). In the ricepile experiment, in addition to the avalanche statistics, another important aspect of SOC systems, i.e., transport properties of sandpiles, was studied through the distribution of residence times (the times spent by the grains inside the pile). Interestingly, the distribution was observed to have a power-law tail and is discussed later.

Simulations
Kertesz and Kiss revisited the BTW's original question of longtime correlations in the toppling activity (avalanche) in the BTW sandpile models by carefully analysing the power-spectrum of the activity time-signal [52]. Kertesz and Kiss showed that the average power spectrum of activity, can be obtained by averaging over the distribution P(s) ∼ s −τ of individual avalanche sizes s. Here it is assumed that an individual avalanche of size s have a Lorentzian power spectrum s 2 /(1 + ω 2 τ 2 s ), with τ s being the mean time-duration of the avalanche of size s; note that we have written the spectrum in terms of the angular frequency ω = 2πf . Provided the power-law dependence of the mean avalanche time duration τ s ∼ s x on avalanche size s, one has S(ω) is the upper cut-off to the avalanche duration. Now, depending on the actual values of the two exponents τ and x, there are two possibilities for the low-frequency behavior of the power spectrum: (i) for 2x + τ > 3, S(ω) ∼ 1/ω (3−τ )/x and (ii) for 2x + τ ≤ 3, S(ω) ∼ 1/ω 2 ∼ 1/f 2 . In simulations, the exponents were found to be τ ≃ 1.1 and x ≃ 0.68, implying that the low-frequency behavior of the power-spectrum is actually 1/f 2 [52]. The theoretical form of the spectrum was verified in direct measurement of the power-spectrum in simulations, which are in fair agreement with 1/f 2 . However, here one should note that the noise spectrum measured in the experiment by Jaeger et al. [28] and that by Kertesz and Kiss [52] (and by BTW in [1]) are in principle two different quantities as the former was related to the outflux at the boundary while the latter to the time-signal of toppling activities in the bulk. In a related work, Jensen et al. [53] pointed out this difference and calculated the power-spectrum of the particle-flow at the boundary. From these studies, one could conclude that the BTW sandpile, though having long-ranged spatio-temporal correlations which are manifest in the powerlaw distributed avalanche size and avalanche duration, apparently does not exhibit the 1/f noise spectrum. The reason because the BTW sandpile does not have 1/f -type spectrum is that the distribution of avalanche duration has a long tail; indeed, if the distribution of avalanche time decays faster, there would have been a possibility of 1/f noise [e.g., provided τ + x = 3, the power spectrum is given by S(f ) ∼ 1/f ]. Manna and Kertesz [54] studied the time series of total mass in the two-dimensional BTW sandpile in the slow driving limit; unlike in the studies in [1,52], here the time unit was chosen as the interval between the addition of two successive grains. The steady-state two-point temporal correlation is the total mass at time t, has an exponential decay, with the correlation time T M scaling with system size L as T M ∼ L 2 . Therefore, the power spectrum of mass fluctuation in the system has a Lorentzian form and essentially has a 1/f 2 decay in the frequency domain f ≫ 1/L 2 . Therefore, the noise activity in the BTW sandpile is less "erratic, " and more correlated, than what was perhaps expected by BTW in their original paper.
In simulations of the BTW sandpiles by Kertesz and Kiss [52], the power spectrum was actually obtained by randomly superposing the activity time-signals during the individual avalanches. As each grain is added only after the system becomes stable (i.e., no unstable site is present), the avalanches created at different positions do not interact, or overlap, with each other. Hwa and Kardar [55] argued that the 1/f noise could arise in the system when the individual avalanches, which are created at different space points, are allowed to overlap with each other and to develop nontrivial correlations in that process. This scenario is possible only when the system is driven at a finite rate, as opposed to the slow driving limit proposed in the theory of SOC. Therefore, according to Haw and Kardar, the BTW sandpile could well exhibit the 1/f noise when the system is driven at a finite rate. Indeed, in many of the observed phenomena, such as voltage fluctuations in conductors, the 1/f noise actually appears at a finite driving rate. If that is the case, presumably one cannot apply the SOC to explain the 1/f spectrum in those situations. Sandpiles driven with a finite addition rate should be applicable to much broader range of phenomena and, moreover, can in principle exhibit interesting long-time correlations. Furthermore, from the physical point of view, another troubling point is that the time unit in the slow driving limit is not quite well defined as the grains are added in unequal time intervals (on the microscopic time scale of the elementary toppling events). Finite driving of the system removes these difficulties and introduces a well defined external time scale in the problem.
Following the update rules of the variants of the BTW sandpiles proposed by Kadanoff et al. [3], Hwa and Kardar considered a one dimensional critical slope-type sandpiles (toppling condition depends on a threshold value of the slope), called "running sandpiles." The system is driven with finite rate where one adds a grain at a site i, i.e., h(i, t + 1) = h(i, t) + 1, with a small probability p = J in /L, where h(i, t) is the height at site i, J in is the average deposition rate and L is the system size. The evolution rules are given by provided that the slope h(i, t) − h(i + 1, t) > (in simulations, n f = 2 and = 8) and with the boundary condition that the left boundary is closed, i.e., h(0, t) = h(1, t) and the right boundary dissipates grains, i.e., h(L + 1, t) = 0. The quantities of interest are the power spectrum of the instantaneous output current J(t) (outflux at the boundary site i = L), which is the number of grains leaving the system at time t, and that of the instantaneous activity or energy dissipation E(t), which is the total number of toppling in the system at time t. In the steady state, the input current balances the average outflux, i.e., J in = J ; note that the model has a maximum output capacity of n f grains per unit time so that J in ≤ n f . The power spectra are found to have the following power-law form S J (ω) ∼ ω −α J and S E (ω) ∼ ω −α E for the outflux and the activity time-signals, respectively. The system is observed to have three distinct scaling regimes, characterized by the different values of the exponents α J and α E : high-frequency single-avalanche regime I, intermediate-frequency interactingavalanche regime II and low-frequency system-wide discharge regime III; each of the regimes are separated from each other by L-dependent time scales [55]. Regime I is dominated by single avalanches and arises in the time range of the typical avalanche duration (as studied in [52], for example). In the simulations, the exponents in the power spectrum are observed to be α J ≈ 2 and α E ≈ 4. Regime II is dominated by the interacting avalanches and happens on the time scale L z , where z is the dynamic exponent. The exponents in the power spectrum in this regime are given by α J ≈ 1 and α E ≈ 1, confirming the presence of the 1/f noise in the system. As we see later, this is the regime (also called hydrodynamic regime), which can be captured by a large-scale hydrodynamic theory [55]. Regime III is dominated by systemwide discharge processes, which happen on the time scales of L 2 , during which the macroscopic slope of the pile changes. Interestingly, the exponents in the power spectrum are negative in this regime and are given by α J ≈ −1 and α E ≈ −1.
The negative exponents suggest anti-correlations in outflux and activity time-signals. It is quite interesting to note that, in a recent study by Garcia-Millan et al. [56], similar anti-correlations, albeit in avalanche time-signals, have been observed in the boundary driven one dimensional Oslo ricepile model [18] in the slow driving limit (grains are added with unit rate and only after the system becomes stable). Here one constructs a time-series of avalanches, where s(t) denotes avalanche size at time step t. Then the variance is observed to be highly suppressed. That is, avalanche time series have hyper-uniform fluctuation, characterized by σ 2 T ∼ T λ with λ = 0. (Normal fluctuations correspond to λ = 1.) In fact, steady-state two-time correlation function Such anti-correlations in time possibly plays a crucial role in determining transport properties of sandpiles and it would be quite interesting to further explore their connections in future.

Theories
At this stage, one would be interested to explore whether some of the above observations can be captured in a unified theoretical framework. Hwa and Kardar [55,57] addressed this issue in developing a large-scale hydrodynamic framework, based on symmetries and conservation laws. Later on, more refined field theories were set up [58], especially in the context of absorbing phase transition with a conserved field, which is associated with the fixed-energy version of sandpiles [20-22, 37, 38, 59, 60]. However, the studies by Hwa and Kardar Frontiers in Physics | www.frontiersin.org made the basis for the standard field theoretical techniques applied to sandpiles, particularly in characterizing the long-time correlations and the corresponding power spectra in various time-signals, such as toppling activities. As we discuss next, the following hydrodynamic theory provides simple predictions of the power spectra of various time-dependent quantities, which could be directly tested in simulations. Motivated by the experiments on granular piles [28] and the observed power spectra of outflux and toppling-activity in simulations of running sandpiles, Hwa and Kardar [55,57] proposed a drift-diffusion equation for height field h(x, t) in sandpiles, which are written in the leading order of non-linearities, where x || = (x.n)n and x ⊥ = x − x || are the position vector along the transport directionn and perpendicular to the transport direction; η(x, t) is a zero-mean non-conserved noise, having , and physically represents random addition of grains; in the conserved case, noise satisfies η( Note that, due to the specific nature of the boundary condition where the grains are added at the one end (say, left) and are dissipated through the other end (say, right) of the pile, the grains tend to flow in a particular directionn. To obtain the above hydrodynamic equation, the local current J(h) is expanded in height variable h and its gradient, where higher order non-linearities, though some of them (e.g., ∇h n and h n , etc. for integer n > 0) permitted by the symmetries, are ignored by assuming the fluctuations are small. It is immediately not clear if the non-linear terms such as h 2 should be allowed as such terms do not obey simple translation in h, i.e., h → h + c. However, Grinstein and Lee [61] later provided a mechanism for generating such terms in the hydrodynamic equations. The absence of this translation symmetry in height can also be seen as the consequence of the boundary condition [e.g., h(x || = L) = 0 at the right end], due to which one cannot translate the height profile arbitrarily. It should be noted here though that the non-linear terms (e.g., that involving λ) are in fact not necessary to generate scale invariance in a system; anisotropic diffusion can be sufficient for that purpose [61][62][63][64].
Provided the stochastic non-linear hydrodynamics Equation ( where J(t) and E(t) are the outflux and the number of toppling at time t. The corresponding power spectra, which are the corresponding Fourier transforms of the time-correlation functions, are given by S J (ω) ∼ ω −α J with α J = 1/z and S E (ω) ∼ ω −α E with α E = 2/z; in two dimensions (z = 6/5), the numerical values of the exponents can be obtained as α J ≃ 0.83 and α E ≃ 1.67, implying the existence of the 1/f -type power spectra in the system corresponding to the anisotropic drift-diffusion Equation (5). For isotropic case (as in the BTW sandpile), one has the simplest non-linear hydrodynamic equation, which, in the absence of noise term η(x, t), is a non-linear diffusion equation, where the local current has the form up to the leading order non-linearities, As mentioned previously, depending on the situation, the noise term can be either nonconservative or conservative. Interestingly, in the above equation for local current, one can see that the bulkdiffusion coefficient (D − λh) can have a clustering instability (the diffusivity reduces for larger value of h). As we see later, the nature of the above diffusive instability is different from that in sandpiles. However, unlike the anisotropic case, the parameters D and Ŵ in the isotropic case can get normalized in a nontrivial way and calculation of the exponents are difficult [58,65]. From an overall point of view, in the dynamic renormalization group analysis, it is not obvious how the renormalized diffusivity could give rise to singular diffusion in the system (the renormalization group theory does predict though "super-diffusive" behavior as z < 2). Hwa and Kardar attempted to provide a large-scale hydrodynamic theory of sandpiles, which could capture scale invariant spatio-temporal structures of the systems having the same symmetries as in sandpiles. The theory was indeed successful in calculating various time-dependent properties of the SOC-like systems, such as the 1/f -type noise spectrum and height correlations, etc. However, the theory remains somewhat unsatisfactory as it cannot explain the precise nature of the transport instability, a unique characteristic of sandpiles, which, as we discuss later, is present in the form of a near-critical singular bulk-diffusion coefficient as well as singular conductivity. Possibly, the large non-linearity due to the threshold-activated dynamics in the SOC systems invalidates the gradient expansion of the local current as given in Equation (6). The current is truncated at the lowest-order non-linearities in height and height gradient, and cannot capture the singular diffusion in sandpiles.

INSTABILITIES IN SANDPILES: SINGULAR TRANSPORT
It is important to characterize instabilities, which can be present in a system near criticality and thus govern the large-scale transport and fluctuations. For example, one could see such an instability in the flow structure of grains in a Hele-Shaw cell where a pile of grains is formed and grown between a confined space and there can be a transition from a continuous flow regime to an intermittent one upon varying certain parameters [66]. In the context of the threshold-activated model systems like sandpiles operating near criticality, the power-law spatiotemporal correlations, as reflected in the avalanche statistics, are strongly indicative of the existence of instabilities in particle transport. However, the precise nature of the transport instability was not clear until Carlson et al. [67] demonstrated the existence of the singular diffusion in a two-state sandpile-like model. Indeed, the bulk-diffusion coefficient in the model was shown to diverge when the system is slowly driven (or, "tuned, " in the sense of fixed-energy sandpiles) toward criticality. But, before we go into this topic, let us first discuss another set of related works, which explored tagged-particle diffusion (also called selfdiffusion) to characterize transport properties of sandpiles.

Tagged-Particle Correlations
In sandpiles, particle transport happens through avalanches. One would therefore expect a close relationship between local toppling activities and transport in the system. As discussed in the previous section, the existence of such a relationship is quite apparent through the emergence of long-time correlations in both outflux and toppling activities (as evident in the respective power spectra), which can be closely related through the height fluctuations in the system (they have similar power spectra with exponents α J ≈ α E ≈ 1). Another way to characterize the timedependent correlations in the system would be to directly connect the toppling activities and the particle transport. This could be done by probing tagged-particle correlations in the systems, e.g., through the distribution of tagged-particle residence times or through characterization of the self-diffusion.
The residence time is defined as the time spent by a tagged particle inside the system. If a particle is added in the system, say, at t in and the particle leaves the system at time t out , the residence time is defined as T = t out − t in . One could ask whether the distribution of residence time is described by a power law and, if so, whether the distribution is any way related to the power-law statistics of avalanches in the system. It turns out that the distribution of residence time, though governed by local toppling activities and can exhibit power-law tails, does not depend on the avalanche exponents. However, the distributions are not universal in the sense that it can depend on the details of particle addition as well as the toppling (or particle-transfer) rules in the systems.
In the Oslo ricepile experiment [29] discussed in the previous section, the distribution of residence times P(T, L) of tracer particles in a system of size L was measured. The tracers were tracked since the time they have been added till the time they leave the system, and the residence times were measured. The distribution was observed to have a power-law tail, P(T, L) ∼ 1/T α , with the exponent α ≃ 2.4 and the mean residence time T ∼ L ν with ν ≃ 1.5. To explain the experimental findings, Boguna and Corral proposed a continuous-time random walk model, having a power-law distribution ψ(t) ∼ 1/t α ′ of trapping time t -the time a grain is trapped or buried in the pile; here the exponent α ′ is an undetermined parameter in the model. The analysis showed that the exponent α in the distribution of residence time is determined by the exponent α ′ in the trappingtime distribution, α = α ′ . Later, using the ricepile model proposed in Christensen et al. [18], it was argued [68] that the power-law tail of the distribution P(T, L) should be dominated by the grains, which are deeply buried in the pile and take a very long time to come out due to the rare height fluctuations. Indeed, the cut-off time T max to the trapping-time distribution is related to the hight fluctuation of the pile and was estimated to be exponentially large T max ∼ exp(κL 3 ) [69], where L is the system size. Moreover, the mean residence time T was exactly shown to be equal to the average mass of the pile, i.e., T ∼ L 2 . The cumulative probability distributions of both trapping times T tr and residence times T have a power-law scaling Prob.(T tr ≥ t) ∼ L ω 1 /[t{ln(t/L ω 1 )} δ 1 ] and Prob.(T ≥ t) ∼ L ω /[t{ln(t/L ω )} δ ], with a logarithmic corrections where ω 1 , δ 1 , ω, and δ are model dependent exponents. Therefore, one obtains α = 2, which also holds in other critical slope-type sandpiles (toppling depends on the threshold value of the slope) as long as one implements "first-in-last-out" rule, i.e., an incoming grain at any site sits on the top, while the topmost grain at the site leaves first [68]. Note that, in the Oslo ricepile experiment, the measured value of α being slightly >2 is possibly because the logarithmic correction was not considered in the experimental data fitting. Though the distribution of residence times in slope-type sandpiles have a universal power-law scaling 1/T 2 (with a non-universal logarithmic correction), the distribution in critical height-type sandpiles (threshold condition is on the height variable) can be a power law or an exponential function [70], depending on the details of the grain addition. Indeed, one can show that the probability distribution (density) function P(x, t) of position x of a tagged particle at time t is governed by a space-dependent diffusion equation, where local diffusivity D(x) of a tagged particle is proportional to n(x), the average number of toppling at position x per unit time; the initial condition can be specified as P(x, t = 0) = r(x) where r(x) is the addition rate at site x [70]. Therefore, unlike in slope-type models where distributions of residence times are determined by the power-law trapping time distributions, the distributions of residence times in heighttype models are governed by the space-dependent diffusivity, which is directly related to the local toppling activity in the systems. We note here that somewhat surprisingly, in both slope-and height-type sandpiles, the exponents characterizing the distributions of residence times are not related to the critical avalanche exponents.
The tagged-particle correlations have also been explored by da Cunha et al. [39,71] in several variants of the conserved Manna sandpiles and both in one and two dimensions, directly through the studies of the self-diffusion coefficient D(ρ) of a tagged particle. Here one considers a stochastic steady-state trajectory of a particle, i.e., position X i (t) of, say, i-th particle at time t. Then the mean-square displacement of the particle up to time t has the following asymptotic form for large t: where the self-diffusion coefficient D(ρ) depends on the number density ρ. Interestingly, it was observed that the self-diffusion coefficient is proportional to the activity in the system. As the activity a(ρ) in the conserved Manna sandpiles has singular behavior near the critical density, the tagged-particle diffusion is also singular in this regime and can be characterized by the same exponent as in the activity (there is no rigorous proof yet though). Indeed, this particular feature is expected to be quite generic in the fixed energy sandpiles undergoing an absorbing phase transition; however, the issue is not settled yet.

Hydrodynamics of Sandpiles
The hydrodynamic equations, and the corresponding dynamic renormalization group (RG) theory, proposed by Hwa and Kardar [55,57] was successful to some extent in explaining the emergence of scale invariance, and in extracting the critical exponents of sandpiles in a general hydrodynamic framework, based on symmetries and conservation laws. Though scope and applicability of this theory is quite broad, the theory cannot really explain the mechanism, which gives rise to the singularity (pole-type) in the diffusion coefficient. The possible reason for the failure to capture the singular diffusion is that the proposed hydrodynamic equations (Equations 5 and 7) do not actually take into account the threshold-activated nature of the local dynamics in the theory. In this scenario, the rigorous derivation of hydrodynamics of the simple two-state model put forward by Carlson et al. [67,72,73] provided not only the much needed impetus into the theory of SOC, but also some useful insights into the possible large-scale structures of sandpiles in general. We discuss below a few examples of simple model systems, where large-scale hydrodynamic time evolution can be obtained and, in some of the systems, the transport can indeed become singular. For simplicity, in the following discussions, we confine ourselves to one dimensional models and conserved (fixed-energy) versions; the models can be suitably generalized to higher dimensions and open boundaries. Such a hydrodynamic framework [74,75] not only provides the insights into the relaxation properties, but can also determine in principle the fluctuation properties of these systems [76][77][78][79].

A Two-State Model of Sandpile
Carlson et al. [67,72] considered a model of lattice gas, consisting of hardcore particles diffusing on a lattice. The occupation variable η i at site i can be at most one: η i = 1 if the site is occupied by a particle, otherwise η i = 0. For simplicity, let us consider the fixed-energy version of the model on a periodic lattice of L sites, where the number of particles is conserved; we denote the number density as ρ = N/L. The model can be appropriately generalized to that with an open boundary and a slow drive. Any particle hops to its nearest vacancy, with a unit rate and symmetrically to its right or the left; the dynamical rules are unlike that in a simple symmetric exclusion process, where a particle can hop only to its vacant nearestneighbor site. For example, consider the following hopping event {. . . 011110 . . . } → {. . . 010111 . . . }, where "1" denotes the hopping particle; note that the reverse process also happens with the same rate. Consequently, the hopping dynamics satisfy detailed balance with respect to Bernoulli product measure, where any site is occupied with probability ρ and is vacant with probability (1 − ρ). Precisely due to the time-reversibility and the product-measure steady state, the model is amenable to a rigorous derivation of hydrodynamics [72], which governs the large-scale relaxation in the system. Here we only present a sketch of the derivation, which will illustrate the procedure of deriving hydrodynamics in a general context. Notably, the model has a long-ranged hopping (avalanche-like) mechanism built into it and, in this way, it mimics the threshold-activated dynamics of sandpiles. Indeed, one can immediately see that, depending on the density ρ, the average jump length l of a particle in one hopping event is given by l = ∞ l=1 lρ l−1 (1 − ρ) = 1/(1 − ρ), which diverges as the density approaches (from below) the critical density ρ c = 1 (i.e., as ρ → 1 − ). The model is not defined for ρ > 1 and one of its variant is discussed in the next section.
To calculate the density-dependent bulk-diffusion coefficient D(ρ), let us consider net local bond-current J(i, i + 1) between i-th and (i + 1)-th sites, Here we have denoted the r-point correlation as A (r) (i) ≡ r r ′ =1 η i+r ′ −r , which is the probability that the consecutive r number of sites to the left of site i are occupied. As one would expect, on hydrodynamic scales, the local density ρ(i, t) = η i (t) is assumed to be a slowly varying function of space and time, and thus the current in the system is small [O(1/L)]. Now, by scaling space x = i/L and expanding A (r) (i+r) ≡ A (r) (ρ(x+r/L)) in the Taylor series around the local density η i (t) = ρ i (t) ≡ ρ(x, t), we have Similarly, the local current can be written using Equation (10), where, in the last step, we have written the correlation A (r) = ρ r by using the product-measure (local) steady-state. Therefore, the rescaled local current J(x) ≡ LJ(i, i + 1) is proportional to the density gradient (rescaled) ∂ρ/∂x and can be written as where the proportionality constant, which is called the bulkdiffusion coefficient, is given by Note that the bulk-diffusion coefficient is in general a non-linear function of density ρ and, in this case, has a third-order pole at density ρ = 1. The time evolution of density at site i is given by can be written, by rescaling space x = i/L and time τ = t/L 2 (called diffusive scaling limit) and in terms of the rescaled or the coarse-grained density field ρ(x, τ ), Provided an initial condition ρ in (x) ≡ ρ(x, τ = 0), the longwavelength density relaxation in the system is exactly described by the above hydrodynamic time-evolution (Equation 15). Carlson et al. [72] also studied tagged-particle correlations, which can provide information about the singularity in the particle transport. Let us denote the position of the i-th tagged particle at time t as X i (t). Then the mean square displacement ( X i (t)) 2 = (X i (t) − X(0)) 2 of the tagged particle has a limiting behavior, where the scaled variance lim t→∞ ( X i (t)) 2 The scaled variance of the displacement of a tagged particle has the order of pole singularity one less than that for the bulk-diffusion coefficient in Equation (14). The bulk and the self diffusivities are in principle two different quantities, but, as discussed above, both of them can be a good indicator of the instability of the particle transport in a system. In fact, we mentioned in the previous section how the tagged-particle diffusion can be used to probe activity in sandpiles. However, it remains to be seen how one can actually connect activity and tagged-particle diffusion on a hydrodynamic level.

A Simple Example
The two-state model discussed in the previous section is not defined above density ρ = 1. Note that, during any hopping move in the model, a particle gets transported to its nearest vacant site instantaneously, introducing possibly an unwanted feature of long-range hopping in the system. Alternatively, to bypass the difficulty, one can introduce a microscopic hopping time-scale into the problem and restrict to only nearest-neighbor hopping [80]. For example, if the number of particles n i (t) at a site i is greater than 1, there is a toppling at unit rate and one particle gets transferred symmetrically to one of its two nearest neighbors. That is, when n i (t) > 1, n i (t + 1) = n i (t) − 1 and n j (t + 1) = n j (t)+1 where j = i±1 is one of the two nearest neighbors of site i. In the modified model, we have relaxed the hardcore constraint on the particle occupancies of the sites. Although this particular model has a trivial scaling behavior, it would illustrate how the relationship between the activity and transport on hydrodynamic scales arise quite naturally in sandpiles.
In fact, given the simple dynamical rules, the model can be treated analytically and the transport coefficients can be calculated exactly as the steady-state, for density ρ > 1, has again a product-measure, with the probability distribution Prob.[n i = n] = P(n) = 1 ρ ρ−1 ρ n−1 of single-site particle-number n = 1, 2, . . . . The above result can be easily obtained by mapping steady-state dynamics of the model (i.e., in the space of steadystate configurations), to that of a zero range process [81]. Indeed, for density ρ > 1, the configurations with n i = 0 are nonrecurrent and does not appear in the steady state; however, in the space of recurrent configurations, the model still satisfies a detailed balance with respect to the steady-state (product) measure. On the other hand, for density ρ ≤ 1, system goes into an absorbing state, devoid of any dynamical activity in the system. One can define an order parameter, called activity, a(ρ) = N a /L, where N a is the number of active sites in the system. It is easy to see that the activity a(ρ) is nonzero for ρ > ρ c = 1 and zero otherwise. For ρ > 1, the activity is obtained exactly as a(ρ) = ∞ n=2 P(n) = 1 − P(n = 1) = (ρ − 1)/ρ, implying a(ρ) ∼ (ρ − ρ c ) β as ρ → ρ + c where order parameter exponent β = 1.
In the presence of a single conserved quantity such as local particle-number, the coarse-grained density field is governed by two transport coefficients-the bulk-diffusion coefficient D(ρ) and the conductivity χ(ρ); both transport coefficients are in general non-linear functions of density ρ. In the context of this simple model, it would be quite instructive to discuss a general theoretical framework to calculate transport coefficients following a linear response theory [77,78,82,83]. The framework provides, in the diffusive scaling limit, an exact hydrodynamic theory of relaxation and fluctuation in the system. As done in the previous section, to calculate the bulk-diffusion coefficient, one has to apply a bias in the form of a small density gradient ∂ρ/∂x and calculate the local diffusive current which is analogous to the Fick's law of diffusion. Likewise, to calculate the conductivity, one has to apply a small external force field of magnitude F, which couples to the local particle-number and thus bias the particle-hopping rates in the direction of the biasing force. The biased or the modified particle-hopping rate c F i→j , from site i to j = i ± 1, is determined using a local detailed balance condition [82,83] as following, where c i→j is the original hopping rate in the absence of the biasing force, m i→j = 1 is the number of particles transferred from site i to j, and δx = 1 is the lattice spacing. The biasing force gives rise to a small drift current which is analogous to the Ohm's law for the electrical conductivity in a current carrying wire. The time-evolution of density ρ i (t) = n i (t) can be straightforwardly obtained using the stochastic update rules for particle number n i (t) in an infinitesimal time interval between t and t + dt as given below, whereâ i is an indicator function for occupancy of ith site (â i = 1 if n i > 0 andâ i = 0 otherwise), δx = 1 is the lattice spacing and The local density evolves according to the equation, which, in the diffusive scaling limit x = i/L, τ = t/L 2 and biasing force F → F/L, leads to the desired hydrodynamic evolution of the scaled density field ρ(x, τ ), where a(ρ) = (ρ − 1)/ρ is the activity, i.e., the probability that a site is active and the local current J(ρ(x, τ )) = J diff + J drift . Now using Equations (16) and (18), we arrive at the expressions of the bulk-diffusion coefficient and the conductivity Note that the bulk-diffusion coefficient and the conductivity both vanish below critical density; however, the diffusivity approaches a constant value while one approaches the critical density from above (ρ → ρ + c ). On the other hand, the conductivity is proportional to the activity in the system and vanishes as density approaches the critical density, and remains zero below critical density. Indeed, as we see next in the conserved Manna sandpiles, vanishing of the conductivity near critical point seems to be a generic feature in conserved stochastic sandpiles, which undergoes an absorbing phase transition. At this stage, one could perhaps recognize the connection to the class of facilitated (or restricted) exclusion processes [84][85][86], where particles hop symmetrically to one of its vacant neighbors, only if the other neighbor is occupied. These facilitated exclusion processes also have similar features in the transport coefficients.
The model discussed above serves as probably the simplest possible example of an absorbing phase transition (APT) [87] in the presence of a conserved quantity. In fact, due to the nonsingular behavior of the order parameter (the order parameter exponent β = 1), the bulk-diffusion coefficient, which is the derivative of the activity w.r.t. density, does not have any divergence. However, as we see below, the situation changes drastically when one modifies the dynamical rules slightly by introducing time-irreversibility (violation of detailed balance) in the steady state, leading to a nontrivial exponent β < 1 and thus singular particle transport in the system.

The Conserved Manna Sandpile
We now consider the conserved Manna sandpile [17], which serves as the paradigm for non-equilibrium absorbing phase transition [23] in a broad class of models with conserved mass, collectively called conserved stochastic sandpiles [32,33,88,89]. Unlike the deterministic bulk dynamics in the BTW sandpile, the update rules in the conserved Manna sandpile, though still governed by a threshold-activated dynamics (critical heighttype), are stochastic in nature. In the conserved version of the Manna sandpile, the total mass in the system remains conserved and the system undergoes an absorbing phase transition below a critical density ρ c . We consider here the continuous-time variant of the Manna sandpile [32], where a site i is selected randomly and, if particle number n i > 1, there is a toppling. During a toppling event, two particles (instead of one particle in the previous model) are transferred stochastically, and independently, to one of its neighboring sites. The hopping move, where two particles go in the opposite directions, is not reversible. As the system undergoes an absorbing phase transition below a critical density ρ c , which in this case is strictly less than 1, the activity has nontrivial behavior as a function of density ρ: while a(ρ) is nonzero for ρ > ρ c , it has a power law decay a(ρ) ∼ (ρ − ρ c ) β , with β < 1, as ρ approaches the critical density from above (ρ → ρ + c ). Clearly, as opposed to the unbounded model considered previously, the activity now develops a singularity as a function of density as the derivative of the activity diverges near criticality. As we see below, the singularity in the activity leads to singular diffusion (as well as singular conductivity) in the system.
To calculate the conductivity, we follow the linear response theory of Chatterjee et al. [78]. We apply a constant external biasing force of magnitude F, which couples to the local particle number and bias the motion in the direction of force according to local detailed balance condition given in Equation (17). In the presence of the biasing force, the stochastic update rules for particle number at a site i can be written as according to Equation (17), the modified hopping rates c F i,0 = 1/2, c F i,+ = (1 + F)/4 and c F i,− = (1 − F)/4 correspond to transfer of one particle to the left and one to the right, that of both particles to the right and that of both particles to the left, respectively. The case with biasing force F = 0 corresponds to the unbiased conserved Manna sandpile. The local density satisfies the following time-evolution equation, Now by scaling space x = i/L, time τ = t/L 2 and biasing force F → F/L, and by assuming a local-steady state [90,91], where activity a i (t) ≡ a[ρ(x, τ )] can be written as a function of local density field ρ i (t) ≡ ρ(x, τ ), one obtains the desired hydrodynamic time-evolution of density field ρ(x, τ ), where local current J(ρ) = −∂a/∂x + a(ρ)F. Comparing with Equations (16) and (18), we get the bulk-diffusion coefficient and the conductivity for the conserved Manna sandpile. One can immediately obtain the near-critical behavior of the bulk-diffusion coefficient using the power-law form of the activity a(ρ) ≃ Const.(ρ − ρ c ) β near criticality, where critical exponent β < 1. Using Equation (28), we indeed get a diverging bulk-diffusion coefficient D(ρ) ∼ 1/(ρ − ρ c ) 1−β → ∞ as ρ → ρ + c [78,92]. On the other hand, the conductivity in the system vanishes as one approaches criticality. Diverging diffusivity and vanishing conductivity near criticality are a unique kind of transport instability, which are usually not seen in equilibrium critical phenomena and could well be the signature of the fixed-energy sandpiles in general.

CONCLUDING REMARKS
We have come a long way since the introduction of selforganized criticality (SOC), and the sandpile model, by Bak, Tang, and Wiesenfeld; it has been an exciting journey, which have opened a new horizon. Indeed, in the past three decades or so, SOC and sandpiles have generated "avalanche-like" activities in the literature, in terms of plethora of exact, experimental and numerical results. Undoubtedly, these studies have helped to shape our understanding of scale invariance in general, and non-equilibrium absorbing phase transition in particular, in a new light. One may however recall the original motivation of BTW, which was to connect the longranged temporal correlations observed in nature in the form of the "1/f " noise to the long-ranged spatial correlations in these systems. In retrospect, it seems that perhaps the BTW's great ambition has not been really fulfilled and the success of achieving the goal of explaining the 1/f noise has been mixed at best. As we have discussed in this article, while, under suitable driving condition, sandpiles do generate 1/f power spectra in time-signal of various transport quantities, in the slow driving limit though it usually generates 1/f 2 -type power spectra.
However, it does not mean in any way that the timedependent correlations in sandpiles are not interesting. On the contrary, sandpiles possess extremely rich spatio-temporal structures, which have some unique aspects in their transport properties and would be worth pursuing for understanding their nontrivial structures better. In the past, the scale-invariant spatial structure in sandpiles have attracted a lot of attention and they have been characterized mainly through avalanche statistics. But, the time-dependent properties are much less explored and the studies in the literature have been perhaps a little scattered. Particularly, our knowledge of the exact hydrodynamic structure of sandpiles concerning long wavelength relaxations is still quite limited, primarily due to the lack of knowledge of the steady (quasi) state measure, which could be partly attributed to the absence of time-reversibility in the systems. Indeed, a comprehensive characterization of the temporal correlations, particularly in terms of the transport properties of the systems, though daunting, would be certainly desirable and could lead to a better theoretical understanding of sandpiles in general.
Moreover we should now emphasize more on the problem of characterizing one of the most interesting aspects-the transport instabilities in sandpiles, which are present in the system near criticality. As demonstrated in the past through the studies of simple model systems, the bulk-diffusion coefficient, which governs the long wavelength density relaxation in the system, diverges near criticality [67,92]. This particular feature, as opposed to that near an equilibrium critical point, where the bulk diffusivity vanishes (known as "critical slowdown"; [93]), is presumably associated with the unique nature of the thresholdactivated dynamics in sandpiles. However, at this stage, it is not quite clear whether the near-critical diverging diffusivity, along with the vanishing conductivity, is a generic feature of sandpiles undergoing an absorbing phase transition and, if so, in what ways the transport instabilities can affect the behavior of the systems near criticality. Understanding of these features can have a far reaching consequence in resolving the long-standing issue of determining the universality classes in sandpiles, which have been hotly debated in recent times [41,60,78,89,[94][95][96][97][98][99][100]. Indeed, apart from the conserved Manna sandpiles, it would be useful to explore the exact hydrodynamic structure in other variants of conserved sandpiles, such as the restricted-height models [88], the conserved threshold-transfer process (CTTP) [34][35][36] or the conserved lattice gases (CLG), etc. [31]; this would help one to check if these models too possess similar instabilities in particle transport and thus to classify a rather diverse set of models accordingly. Moreover, an exact hydrodyanamic characterization of sandpiles can provide insights not only into the relaxation properties, but also the fluctuation properties of the systems [78].