{'fraction_non_alphanumeric': 0.03944423722491608, 'fraction_numerical': 0.008718761656098471, 'mean_word_length': 4.213660670879922, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Continuously varying exponents in a sandpile model with dissipation near surface

We consider the directed Abelian sandpile model in the presence of sink sites whose density f_t at depth t below the top surface varies as c~1/t^chi. For chi>1 the disorder is irrelevant. For chi<1, it is relevant and the model is no longer critical for any nonzero c. For chi=1 the exponents of the avalanche distributions depend continuously on the amplitude c of the disorder. We calculate this dependence exactly, and verify the results with simulations.


Introduction
An important question in the general area of self-organized criticality [1] is the universality of critical exponents. Do small perturbations of the local evolution rules leave the critical exponents unchanged? How do perturbations of the system at the boundary affect the critical steady state? These are the questions we study in this paper, and construct a simple analytically tractable model where changes near the boundary lead to continuously varying critical exponents.
The influence of boundary perturbation in nonequilibrium systems is more complicated than in equilibrium systems [2]. In the latter, away from phase transition points, the bulk is largely unaffected by the boundary, except that in case of spontaneous symmetry breaking, the boundary conditions will remove the degeneracy between different phases, and some finite-size effects. Near critical phase transitions, the critical behavior near the boundary is substantially modified, and constitutes a different class of (surface) critical phenomena, characterized by a new set of (surface) exponents [3]. In the so-called 'ordinary transitions', the surface phase-transition is driven by the bulk transition in the sense that the transition occurs at the critical temperature corresponding to the bulk, and its position or the critical exponents are not affected by local changes in the Hamiltonian near the surface.
In self-organized critical systems, the correlation length is, by definition, always large, and the influence of perturbations at the boundary is felt deep inside. In the well-known Bak-Tang-Wiesenfeld (BTW) sandpile model of self-organized criticality [4], even the existence of a critical steady state depends on the presence of a surface where particles can leave the system. Thus, in such systems, boundary perturbations affect the overall behavior of the system to a much greater extent. As the steady state is critical, various measured average quantities in the steady state near a surface will differ from the bulk values by an amount which will decrease as a power law of the distance from the surface.
For example, in the BTW model in two dimensions, Brankov et al. [5] have shown that the mean height of the pile at a site at distance t from the surface is less than the mean value in the bulk by an amount proportional to t −2 . In a one dimensional sandpile model with threshold dissipation, Ali found that the density is larger near the boundary [6] and the excess density varies as 1/t. We would like to understand the effect of such extended perturbations on other properties of the system. For example, one may ask if such a density profile significantly affects the probability of extinction of avalanches that reach the boundary. One could even ask if it determines the avalanche exponents [7]. In this context, it seems worthwhile to understand the properties of self-organized critical systems with 'artificially produced' power-law surface perturbations.
In this paper, we consider the directed Abelian sandpile model in the presence of some sink sites, whose density varies with the distance from the surface with a power χ. We are able to determine the χ-dependence of the avalanche exponents exactly. Our analysis shows that the disorder is irrelevant if χ > 1, and is relevant for χ < 1. For χ = 1, the perturbation is marginal, and the critical exponents depend continuously on the disorder strength. In the case of equilibrium critical phenomena, for Ising model with deviation of the coupling strength from the critical value varying with distance from surface, a similar behavior is seen [8], and may be understood in terms of a simple renormalization argument [9].
In fact precisely this question comes up in the sandpile model studied by Tadić and Dhar (TD) [10]. This model (its definition is given later in the paper) belongs to the universality class of directed models with stochastic toppling rule, with no multiple topplings [11]. For this model, TD presented some intuitive arguments to determine the exact critical exponents for a general dimension d in terms of the exponents of the directed percolation problem. The arguments, involving determining the statistics of clusters in a directed percolation model with space-varying density, are plausible, but not rigorous. In the model studied here, our exact answer coincides with that obtained based on the same reasoning, thus provides some additional evidence in the latter's favor.
In the next section we define our model. We determine the disorder-averaged two-point correlation function of toppling events exactly in section 3. This shows the qualitative difference in the behaviors of the avalanche distribution for χ > 1 and χ < 1. We give heuristic arguments to determine the avalanche exponents as a function of the amplitude c in the marginal case χ = 1. These results are confirmed by numerical investigations. A somewhat more careful derivation of the avalanche exponents requires the consideration of the three-point correlation function of toppling events which is presented in section 4. The connection to the argument used in [10], and some generalizations to other systems are discussed in the concluding section.

Definition of the model
We consider a two dimensional lattice in the shape of an open cylinder (see Fig. 1). Each site is labelled by two integer coordinates (x, t), where x labels the coordinate in the horizontal direction, and t denotes the depth of the layer (0 ≤ t ≤ T ). It is convenient to work with only the set of sites for which x + t is even (one of the two sublattices into which this lattice may be decomposed). We denote the number of sites in each horizontal row by L. To each lattice site (x, t), a non-negative integer variable h(x, t) is associated, called the height of the sandpile at that site. A sandgrain is added at a randomly chosen site x on the top layer (t = 0). If the height of a site (x, t) exceeds the critical value h c = 2, i.e, if h(x, t) ≥ 2, it is unstable and it topples. On toppling at a site, its height is decreased by 2, and one grain is transferred to each of the two downward neighbors (x ± 1, t + 1). In this way the transferred sand-grains may activate the two downward sites and an avalanche of relaxation events may take place. Some typical avalanches are shown in Fig. 2.
The size of an avalanche may be measured in terms by its downward length (called its duration t), or the number of toppled sites (called its mass m). In the critical steady state, the probability distribution of avalanche duration and avalanche mass has power-law tails [12], i.e., In this paper, we consider the effect of introducing surface disorder in the model. We do this by declaring a small fraction of sites as sink sites. Any particle falling onto a sink site is lost from the system. The sink sites are uncorrelated with each other. We assume that the concentration of sink sites f t at the depth t decreases as a power of t: The case where the sink concentration is uniform was studied earlier by Tadić et al. [13]. They showed numerically that any finite concentration of sinks destroys the criticality of the steady state. Theiler [14] showed that in the mean-field theory a uniform sink density f induces an exponential cutoff in the depth reached by avalanches which varies as f −1 for small f .
In the next sections we investigate the two-point and three-point correlation functions which allow to derive the avalanche exponents exactly.

The two-point correlation function
The characterization of the steady state is very simple in our model. From the general theory of Abelian sandpile models [15], it is easy to see that all stable configurations are recurrent, and occur with equal probability in the steady state. Thus, in the steady state, each non-sink site is empty or occupied with equal probability, independent of the occupation of other sites.
Let us now consider the propagation of avalanches in the system. Let G(x, t|x ′ , t ′ ) be the probability in the steady state that a toppling will occur at (x, t) if a particle is added at (x ′ , t ′ ). If the site (x, t) is a sink site, we define G(x, t|x ′ , t ′ ) to be zero. In the steady state, at non-sink sites, the average number of particles coming into (x, t) must equal the number of leaving particles. Thus, the two-point correlation function G(x, t|x ′ t ′ ) satisfies the equation where η(x, t + 1) is zero if the site is a sink site, and one otherwise. Averaging over the configuration of sink sites gives the disorder-averaged two-point correlation for all x 0 . We use Eq.(5) as a recursion equation and determine G(x, t|0, 0) for all t, given its initial value at t = 0 It is easy to show that the disorder-average function G is given by where G 0 is the propagator for the problem with no sinks, i.e., and where P t is defined to be the product For the case with no disorder, it is known that Prob{duration ≥ t} and G 0 (0, t|0, 0) both vary with the same power of t [12]. If we make the plausible assumption that this continues to hold in the presence of sinks, we expect that the probability distribution of the avalanche duration is of the form since the propagator G 0 is known to decrease asymptotically as t −1/2 for large t [12]. Let us consider some limiting cases. In the case of homogeneous disorder f t = f we get that P t ∼ (1 − f ) t , and the probability that an avalanche reaches depth t decreases exponentially with t. This is in agreement with the numerical results of [13].
In the case that the disorder density f t decreases as a power-law according to Eq. (3) one has to distinguish three cases: i) For χ > 1, the disorder is irrelevant. P t tends to a finite constant value for t → ∞, and the asymptotic behavior of the disorder averaged correlation function G is the same as that of G 0 , i.e., the disorder does not affect the scaling behavior. A typical avalanche for χ = 2 is shown in Fig. 2. While branchings can occur, most of the longer avalanches have only one prominent branch.
ii) For χ < 1, the disorder is relevant, and leads to loss of criticality in the system. P t tends to zero as exp (−Kt 1−χ ), where K is a constant which depends on χ and c. Thus the function G decreases to zero with the distance as a stretched exponential. In contrast to the case χ > 1 the avalanche propagation for χ < 1 is characterized by many branchings, though only one branch survives for long (see Fig. 2).
iii) For χ = 1, the disorder is marginal and the avalanche distribution exponent varies continuously with the parameter c. Here, P t tends to zero for large t as t −c . From Eq. (8), we see that G decreases as t −1/2−c and thus we find The avalanche distribution Prob{duration ≥ t} still exhibits a power-law decay, but the exponent depends continuously on the amplitude c. This behavior has been verified in our simulations. Figure 3 displays the probability distribution Prob{duration ≥ t} for three different values of c. For sufficiently large durations the asymptotic behavior agrees with Eq. (12). Let us now consider the distribution of masses of avalanches. We see from the explicit form of the avalanche propagator G that while its dependence on t is modified by the disorder, its dependence on the transverse coordinate x is unaffected. Hence, even in the presence of disorder, the average transverse size of an avalanche at depth t scales as t 1/2 . As each site in an avalanche topples only once, the average width of an avalanche cluster varies as t 1/2 and the total number of toppled sites m scales as t 3/2 . Then Prob{mass ≥ m} ∼ m −τm+1 with In our simulations, we used lattices of size L × T with T = L. In this case, the transverse size is effectively infinite as even the largest avalanches are not affected by the finite transverse size. Thus the finite-size effects are governed by the longitudinal size (T ) alone, and the mass distribution is expected to fulfill the finite-size scaling ansatz [16] Prob If finite-size scaling works a data collapse of the distributions for different system sizes L×T (with L = T ) is obtained for ν = 3/2 and β = 1/2+c. The corresponding scaling plots for three different values of c are shown in Fig. 4. In all cases we observed good data collapse which confirms the above analysis.

The three-point correlation function
In the previous section we obtained the avalanche exponents by using some scaling arguments, which are not careful in distinguishing the probability that a given site at depth t is toppled in an avalanche from the probability that at least one site at depth t is toppled. We now present a cleaner derivation of the above results.
Let us denote the number of sites toppling in the layer t by N (t). Based on the experience with the model without disorder, we make the assumption that the distribution of avalanches in this model can be described by the scaling function for n = 0 and where the scaling function F(x) decreases fast for large x, and A is a constant. From this scaling form, it follows that N (t) ∼ t 2b−a , and N 2 (t) ∼ t 3b−a . Thus knowing the t-dependence of N (t) and N 2 (t) , we can determine the exponents a and b. Clearly, using Eq.(8) we get We would now like to calculate N 2 (t) . Let G (3) (x 1 , x 2 , t) be the probability that both sites x 1 and x 2 at depth t topple when a particle is added at x = 0 on the top layer t = 0, averaged over all configurations of the disorder. Then, clearly, we have For t > 0, the three-point correlation function G (3) satisfies the equation for x 1 = x 2 . We have to supplement this equation with the boundary condition for As for the case with no disorder, Eq.(18) is solved by the ansatz which satisfies Eq.(18) automatically for x 1 = x 2 . We choose the function f (y, u) so that the equation holds also if x 1 = x 2 . Putting x 1 = x 2 = x in the previous equation, we get This is a linear equation connecting f (y, u) to its values at sites higher up. Hence these can be solved recursively. The calculation of N 2 (t) can be further simplified. We define the function and note that where we define P −1 = 1. Then using Eqs. (17,20), we get The recursion equation for F (u) is also easily deduced from the recursion equation for the function f (y, u). We note that x G 2 (x, t|y, u) is independent of y and only depends on t and u. Let us write Then summing Eq.(21) over x, we get The effect of the disorder is just to introduce an overall rescaling factor to K(t, u) in the absence of disorder: This implies that the asymptotic behavior of F is given by Substituting in Eq.(24), we see that The obtained scaling behavior of N (t) and N 2 (t) implies that the exponents of the assumed scaling form [Eq. (15)] are given by a = 1 + c, and b = 1/2. This justifies the earlier heuristic calculation of the exponents. It is straight-forward to extend this treatment to higher dimensions. We only mention the result here: in two (i.e. 2 + 1 dimensions) and higher dimensions the avalanche exponents are given by τ t = 2 + c, and τ m = (3 + c)/2.
As in the case with no disorder, there are logarithmic correction terms to the exponent values given above in the two-dimensional case.

Discussion
Let us comment briefly on the relationship of this model to that studied earlier by TD. Their stochastic variant of the directed sandpile model is defined on the same lattice as presented in Fig. 1, and particles are added on the top layer at random, and removed from below. The difference from the model studied here is that there is no quenched randomness, and the stochastic toppling rules are defined as follows: If a particle is added to a site, only then is the site examined for toppling. If the number of particles exceeds 2, with a non-1 probability p a toppling occurs and one particle is transferred to each of the two downward sites, and with probability (1 − p) nothing happens. While there is particle conservation in the TD model, so far as the evolution of a single avalanche is considered, the sites where no toppling occurs when particles are added, act as sink sites. The density of sink sites is higher near the top in both models. The randomness in the TD model provided by the stochastic toppling rules is 'annealed', unlike the quenched randomness in the model studied here.
In [10], it was shown that this model has a steady state for p ≥ p ⋆ , where p ⋆ is the critical probability for directed site percolation on a square lattice. In the steady state, there is a small fraction of sites having number of particles 2, 3, 4, . . ., but the average density is finite, and is a function of the depth of the layer.
TD argued that a density less than the bulk density locally corresponds to a finite correlation length which is spatially varying. They argued if the medium has an inverse correlation length κ(t), which is small, and varies slowly as a function of t, then the survival probability of a disturbance up to depth t is given by The first factor is what we should get if κ(t) is zero everywhere. The second term is a simple generalization of the exp (−κt) factor when κ(t) varies with t. We note that the expression is non-perturbative, in that the small parameter κ is exponentiated. More importantly, κ is usually a nontrivial power of the starting perturbation δρ(t) in the density (in the TD case, κ(t) ∼ [δρ(t)] ν , where ν is the longitudinal correlation length of directed percolation problem). In the TD model, κ(t ′ ) varies as c/t ′ , and hence the argument of the exponential varies as −c log t, for large t. This immediately gives the c-dependence of the power-law behavior of Prob(t).
Using explicit calculation of correlation functions for small finite t, it is easy to verify that Eq. (31) is not exactly valid for the TD case. However, for the subsequent derivation the exponents of the TD model to be valid, the ansatz has to be at least asymptotically exact. As a proper treatment of the correlations in the steady state of TD model seems very difficult, one may try to find other simpler models where the validity of the ansatz can be demonstrated.
The steady state in our model is very simple, and thus the model is a good testing-ground. It is gratifying that, in the present model, the ansatz does become exact, and the two-point correlation function satisfies Eq. (8) exactly. Thus we have shown that this approximation is exact for our model.
The approximation is in the spirit of the familiar semiclassical Wentzel-Kramers-Brillouin-Jeffreys approximation in quantum mechanics [17]. We note that the mathematical mechanism giving rise to continuous variation of critical exponents with amplitude of perturbation applies to a much wider class of problems. Consider, for example, the generalization of the equilibrium problem treated in [8] involving extended surface defects to general p-vector spins {S(n)} in ddimensional half-space. Lattice sites are labelled by a d-dimensional integer vector n = (n 1 , n 2 , . . . , n d ) with n 1 > 0. The Hamiltonian of the system is H = J(n, n ′ ) S(n) · S(n ′ ) where the summation is over all nearest neighbor pairs of sites n, n ′ . The coupling constant J(n, n ′ ) is assumed to depend only on n 1 (we may assume that n 1 ≤ n ′ 1 ) as where J ⋆ is the value of coupling in the bulk, and C > 0.
As mentioned earlier, the simple renormalization argument due to Burkhardt [9] shows that for χ > 1/ν, where ν is the (bulk) correlation length exponent, this perturbation is irrelevant, and for χ < 1/ν, it is relevant. In the marginal case χ = 1/ν, we get the inverse correlation length κ(n 1 ) ∼ C ν /n 1 . Let us consider the spin-spin correlation function S(n) · S(n ′ ) . For simplicity, we may take n and n ′ to differ only in the first coordinate. Then, it is easy to write an expression for this correlation function using the ansatz Eq. (31). The mathematical problem then is the same as discussed in this paper, and one finds that in the marginal case χ = 1/ν, and the critical exponent corresponding to the decay of this correlation function varies linearly with C ν .
Of course, the dramatic effect of the boundary-perturbations on the exponents of the avalanche size distribution is due to the fact that all avalanches are initiated at the top boundary in our model. If we add particles at a fixed distance t 0 from the boundary, it is easy to see that large avalanches of duration much greater than t 0 have the same asymptotic behavior as the case t 0 = 0. However, shorter avalanches find significantly fewer defects than if they had been added at the top. If we added particles everywhere in the bulk, the avalanche size distribution is averaged over different values of t 0 , and in the thermodynamical limit is the same as the case without defects.
We thank M. Barma for a critical reading of the manuscript. S. Lübeck would like to thank V. B. Priezzhev and A. Hucht for useful discussions.