The abelian sandpile model on a random binary tree

We study the abelian sandpile model on a random binary tree. Using a transfer matrix approach introduced by Dhar&Majumdar, we prove exponential decay of correlations, and in a small supercritical region (i.e., where the branching process survives with positive probability) exponential decay of avalanche sizes. This shows a phase transition phenomenon between exponential decay and power law decay of avalanche sizes. Our main technical tools are: (1) A recursion for the ratio between the numbers of weakly and strongly allowed configurations which is proved to have a well-defined stochastic solution; (2) quenched and annealed estimates of the eigenvalues of a product of $n$ random transfer matrices.


Introduction
The abelian sandpile model (ASM) is a thoroughly studied model both in the physics and in the mathematics literature see e.g. [3,12,9,15] for recent review papers on the subject. In physics, it serves as a paradigmatic model of self-organized criticality (SOC). SOC is usually referred to as the phenomenon that the model exhibits power law decay of correlations or avalanche sizes, without fine-tuning any external parameters such as temperature or magnetic field. In mathematics, the ASM is connected to several combinatorial objects such as spanning trees, graph-orientations, dimers, and it has an interesting abelian group structure.
The ASM has been studied on the Bethe lattice (i.e., the rootless binary tree) in [5]. Via a recursive analysis, based on a transfer matrix method, the authors in [5] arrive at exact expressions of various quantities of interest, such as the single height distribution, correlation functions of height variables, and avalanche size distribution.
There are various motivations to consider the ASM on random graphs. As an example, we mention integrate-and-fire models in neuroscience (see e.g. [10] and related papers), where the connections between neurons are updated after a neuron has fired. The typical connection structure of a network of firing neurons is therefore generically not translation invariant, and time dependent. As a first approximation, one can quench the randomness of the connection graph and study the firing of neurons on the derived random graph. So far, the ASM model has been studied on small world graphs from the physicist's perspective, using a renormalization group approach [8]. In the mathematics literature, there are recent studies on so-called "cactus" graphs [13].
In this paper, we start this study of the ASM on random graphs with the ASM on a random tree, for the sake of simplicity chosen to be a realization of a binary branching process with branching probability p ∈ [0, 1]. We use the transfer-matrix method of [5] to express relevant quantities such as the correlation of height variables and the avalanche size distribution in terms of the eigenvalues of an ad hoc product of random matrices. This is the fundamental difference between the Bethe lattice case and the random tree, namely the fact that the transfer matrices depend randomly on the vertices and instead of having to deal with the n-th power of a simple two by two matrix, one has to control the product of n random matrices.
The crucial quantity entering the transfer matrices is the so-called characteristic ratio, which is the ratio between the numbers of weakly and strongly allowed configurations. This ratio is equal to 1 in the infinite Bethe lattice for every vertex and it is close to 1 for vertices belonging to a finite subset of the Bethe lattice which are far away from the "boundary" (see later on for precise statements). In our case, we show that for an infinite random tree the characteristic ratio is a well-defined random variable, uniquely determined by a stochastic recursion. The transfer matrices will then contain elements with that distribution. We also consider deterministic trees that are strict subsets of the binary tree where the characteristic ratio can be computed explicitly. Next, we prove the exponential decay of correlation of height variables (as in the Bethe lattice case), and show that for a branching probability p sufficiently small, but still supercritical (p > 1/2), i.e., the branching process survives with positive probability, avalanche sizes decay exponentially. This shows a transition between exponential decay of avalanche sizes, for p small, and power law decay for p close to (possibly only equal to) one.
Our paper is organized as follows. First we recall some basic material about the ASM on trees and the recursive technique developed in [5]. Second we study the recursion for the characteristic ratio and show it has a unique solution for the random binary tree. Finally we give quenched and annealed estimates of the eigenvalues of the product of n random matrices, which we apply in the study of correlation of height variables and avalanche sizes.

Abelian sandpile model on subtrees of the full binary tree
We summarize here the basic and standard objects of the abelian sandpile model on a (general) tree. More details can be found e.g. in [5,11].

Rooted and unrooted random trees
We denote by B n the rooted binary tree of n generations, and by B ∞ the rooted infinite binary tree. For a more general tree T we write T i if we want to indicate that the tree has root i. The rootless infinite binary tree or Bethe lattice is then obtained by joining two infinite rooted binary trees by a single edge connecting their roots. A random binary tree of N generations with branching probability p ∈ [0, 1] is a random subset T N of B N obtained as follows. Starting from the root, we add two new vertices, each connected with a single edge to the root (resp. no vertices), with probability p (resp. (1 − p)), and we iterate this from every new vertex independently for N generations. By letting N → ∞ we obtain the full binary branching process. Joining two independent infinite copies of this process by a single edge connecting their roots creates the rootless random binary tree. This last procedure is of course identical to create the non-random rootless binary tree from non-random rooted binary trees.

Height configurations and legal topplings
For T a finite subtree of the Bethe lattice, height configurations on T are elements η ∈ {1, 2, . . .} T := H T . For η ∈ H T and u ∈ T , η u denotes the height at vertex u.
For a configuration η ∈ H T , we define the toppling operator T u via where ∆ is the toppling matrix, indexed by vertices u, v ∈ T and defined by (u, v neighbors in T means that an edge of T connects u to v). In words, in a toppling at u, 3 grains are removed from u, and every neighbor of u receives one grain. A toppling at u ∈ T is called legal if η u > 3. A sequence of legal topplings is a composition T un •. . .•T u 1 (η) such that for all k = 1, . . . , n the toppling at u k is legal in T u k−1 • . . . • T u 1 (η). The stabilization of a configuration η ∈ H T is defined as the unique stable configuration S (η) ∈ Ω T that arises from η by a sequence of legal topplings.

Addition operator and Markovian dynamics
For T a finite subtree of the Bethe lattice, and for u ∈ T , the addition operator is the map a u : Ω T → Ω T defined via where δ u ∈ {0, 1} T is such that δ u (u) = 1 and δ u (z) = 0 for z ∈ T , z = u.
In other words, a u η is the effect of an addition of a single grain at u in η, followed by stabilization. The addition operators commute, i.e., a u a v = a v a u . This is the wellknown and crucial abelian property of the sandpile model.
The dynamics of the sandpile model is then the discrete-time Markov chain {η(n), n ∈ N} on Ω T defined via where X i are i.i.d. uniformly chosen vertices of T . Given a stable height configuration η and u ∈ T , we define the avalanche Av(u, η) induced by addition at u in η to be the set of vertices in T that have to be toppled in the course of the stabilization of η + δ u .

Recurrent configurations and stationary measure
The recurrent configurations of the sandpile model form a subset of the stable configurations defined as follows. A configuration η ∈ Ω T contains a forbidden subconfiguration (FSC) if there exists a subset S ⊂ T such that for all u ∈ S, the height in u is less than or equal to the number of neighbors of u in S. The restriction of η to S is then called a FSC. A configuration is allowed if and only if it does not contain a FSC. Recurrent configurations coincide with allowed ones, and are collected in the set R T .
We denote by P(Ω T ) the set of probability measures on Ω T . The Markov chain (3) has a unique stationary probability measure µ T ∈ P(Ω T ) which is the uniform measure on the set R T where δ η is the point mass concentrated on the configuration η.

Specific properties of the sandpile model on a tree
In this section, for the sake of self-containdness, we briefly summarize some basic facts from the paper [5] which we need later on. In a subset of the Bethe lattice, the distance between two vertices is defined as the length (i.e., the number of edges) of the shortest path joining them. A vertex is a surface vertex if it has a number of neighbors strictly less than 3.

Weakly and strongly allowed subconfigurations
The class of allowed configurations can be divided into weakly and strongly allowed ones. Let T be a rooted finite tree with root u and extend it with one vertex v and one edge < u, v >. Consider an allowed configuration ξ on T . We put height 1 at vertex v and investigate the derived configuration ξ on T ∪ {v} such that ξ is the restriction of ξ to T , and ξ v = 1. If ξ has no FSC, then we call ξ strongly allowed on T , otherwise weakly allowed. We give in Figure  1 an example of weakly and strongly allowed configurations. On one hand, if η u = 2, then there exists a forbidden subconfiguration on S = {v, u, i 1 }.
On the other hand, if η u = 3, there are no forbidden subconfigurations.

Characteristic ratio and recursion
Let T be a finite tree, rooted or not. A key quantity in the analysis of [5] is the characteristic ratio x(T ) ∈ [1/2, 1] of T between weakly and strongly allowed configurations. For the empty tree T = ∅ we put x(∅) = 0.
For an infinite tree T we say that the characteristic ratio is well-defined if the limit lim T ↑T x(T ) exists, where lim T ↑T is taken along the net of finite subtrees. That is, lim T ↑T x(T ) = a means that for every ε > 0, there exists The characteristic ratio satisfies a recursion property both for rooted and unrooted trees: where T 1 , T 2 are two non-intersecting subtrees of the tree T defined as follows.
If T = T o is a rooted tree with root o, then T 1 , T 2 are obtained by deleting the root o and splitting the tree into two subtrees whose roots are the two descendants of the root o, see Figure 2. If T is an unrooted tree, we pick one Figure 2: Example of a tree with root o which splits into T 1 and T 2 of its edges, < u, v >, we delete this edge and split T into two rooted subtrees T 1 and T 2 whose roots are u and v, see Figure 3. This recursion (4) holds for Figure 3: Example of a rootless tree which splits into T 1 and T 2 finite trees, and by passing to the limit for infinite trees, provided the limits defining x(T ), x(T 1 ), x(T 2 ) exist. Dhar & Majumdar obtain from (4) that when the distance from the root of T to the nearest surface vertex tends to infinity (what they call "deep in the lattice"), the characteristic ratio x(T ) tends to 1. The authors then derive several explicit quantities (such as height probabilities of a vertex u deep in the lattice) by replacing the characteristic ratio by 1.

Transfer matrix approach
The transfer matrix approach allows to compute the two-point correlation functions and avalanche size distribution. Let u, v be two vertices in the tree at mutual distance n. They determine n + 3 subtrees T 1 , . . . T n+3 , see Figure  4. The number of allowed configurations when fixing the heights at u and v can then be obtained via a product M of n + 3 two by two matrices with elements determined by the characteristic ratios of the subtrees T 1 , . . . , T n+3 (see Section 4 below for the precise form of these matrices). In particular ... Figure 4: Example of a tree with n + 3 subtrees for the Bethe lattice, all the matrices involved in this product are equal to 2 2 1 3 because the characteristic ratios of all the involved subtrees are equal to one. Let T be a finite or infinite subtree of the full binary tree. The two-point correlation function, i.e., the probability that two vertices u, v at mutual distance n have height i resp. j, for i, j ∈ {1, 2, 3}, is equal to (see [5, Section 5, eq. (5.11)]) where λ − (M )/λ + (M ) is the ratio of the smallest and largest eigenvalues of the matrix M and a i,j are some numerical constants depending on i, j. For the binary tree, λ − (M ) = 1 and λ + (M ) = 4 n . The avalanche size distribution is determined by the inverse of λ + (M ). For the binary tree, upon addition of a grain at a vertex u, the probability that the avalanche Av(u, η) is a given connected subset C of cardinality n containing u is equal to for some constant C, independent of the shape of the subset C . Since there are 4 n n −3/2 (1 + o(1)) connected subsets of the Bethe lattice of cardinality n containing u, one concludes that for large n (see [5, Section 6, eq. (6.13), (6.14)]), i.e., the tail of the avalanche size distribution decays like n −3/2 .

Some characteristic ratios
On the Bethe lattice x(T i ) = 1 for the (infinite) rooted subtrees (i = 1, 2, 3) attached to every vertex since every T i is an infinite rooted binary tree (see Subsection 2.5.2). This property considerably simplifies the analysis of [5] and is no longer valid in the inhomogeneous or random cases studied here.

The characteristic ratio of the random binary tree
Let T n be the random binary tree of n generations starting from a single individual (the "zero-th" generation) at time n = 0 (cf. Section 2.1). Then for n > 0, x(T n ) satisfies the recursive identity [5, Section 3, eq. (3.12)] and T i n−1 are the (possibly empty) subtrees emerging from the (possibly absent) individuals of the first generation. For a tree T 0 consisting of a single point we have because heights 2, 3 are strongly allowed and height 1 is weakly allowed (cf. Section 2.5.1). This value 1/2 can also be obtained from the recursion (8) by viewing a single point as connected to two empty subtrees (for which and analogously in the other argument.
PROOF. The proof is straightforward and left to the reader.
PROOF. Denote by µ n ∈ P(Ω Tn ) the distribution of x(T n ). The recursion (8) induces a corresponding recursion on the distributions µ n+1 = F (µ n ). We show that F is a contraction on the set P([1/2, 1]) of probability measures on [1/2, 1] endowed with the Wasserstein distance. This implies that it has a unique fixed point µ * and from every initial µ 0 , µ n → µ * in Wasserstein distance and thus weakly. Let g be a Lipschitz function. Then we have where T 0 , T 1 , T 2 are the three trees which cannot be split into two subtrees both non-empty, see Figure 5. Every other tree T (appearing in the last sum of (12)) can be split into two subtrees both not reduced to a single point. Then using (8), the expression in (12) becomes is a contraction on P([1/2, 1]) endowed with the Wasserstein distance. The contraction factor is bounded from above by 8/9.
PROOF. Denote by L the set of Lipschitz functions g : [1/2, 1] → R with Lipschitz constant less than or equal to one, i.e., such that |g(x) − g(y)| ≤ |x − y| for all x, y. We use the following two formulas for the Wasserstein distance of two elements µ, ν of P([1/2, 1]) [6]: where in the last right hand site the infimum is over all couplings P with first marginal P 1 (resp. second marginal P 2 ) equal to µ (resp. ν).
To estimate d(F (µ), F (ν)) we start with the first formula (15) Now use the definition of f and the fact that x, y ∈ [1/2, 1] to estimate This gives, using the Lipschitz property of g and a coupling P of µ and ν: Taking now the infimum over all couplings P, using (16) to bound (18), and taking the supremum over g ∈ L yields sup g∈L A(µ, ν, g) ≤ 1 4 d(µ, ν).
To estimate the term B(µ, ν, g), use the elementary bound (since x, x , y, y ∈ [1/2, 1]): We then have, using the Lipschitz property of g, and a coupling P = P (dxdx ; dydy ) of µ ⊗ µ and ν ⊗ ν Taking the supremum over g ∈ L , and infimum over the couplings P , we find sup Combining (19), (20) with (17) we arrive at where in the final inequality we used the elementary bound This proof of Lemma 3.2 completes the proof of Proposition 3.1.
REMARK 3.1. The random binomial tree is such that every vertex has two children with probability p 2 , 1 child with probability 2p(1 − p) and 0 children with probability (1 − p) 2 , p ∈ [0, 1]. For this case the same idea as in Lemma 3.2 gives a recursion leading to a contraction in the Wasserstein distance and hence x(T n ) has also a unique limiting distribution.

The characteristic ratio of some deterministic trees
The recursion (8) allows also to compute the characteristic ratio for certain (deterministic) infinite subsets of the full binary tree in terms of iterations of sections of f .

Infinite branch
First, consider a single branch of length n, i.e., the tree T min n consisting of a root and n ≥ 1 generations of two individuals each (see Figure 6). Using is monotonically increasing in n with limit which is the unique positive solution of REMARK 3.2. We can generalize this "backbone" tree to a general "backbonelike" tree where at each point the same finite tree T is attached (it is a singleton in the backbone case). The characteristic ratio equals the positive solution of the fixed-point equation where x(T ) denotes the characteristic ratio of the tree T , i.e., x * = 1 2 (−1 + 5 + 4x(T )).

Finite perturbations of a single branch
Attaching a finite subtree T of the full binary tree B ∞ at level n in the infinite single branch T min ∞ leads to a tree T per n with characteristic ratio where ϕ(x) = f (1/2, x) (cf. (23)) is applied n times. This shows that inserting a finite tree T at level n has an effect on the characteristic ratio that vanishes in the limit n → ∞, exponentially fast in n. Moreover, since x(T ) ≥ 1/2 and x → f (x, y) is monotone for all y (by Lemma 3.1), we have from (26) From (27) and (22) we conclude that for every infinite subtree T ∞ ⊂ B ∞ for which x(T ∞ ) exists,

Transfer matrix and eigenvalues: uniform estimates
In the analysis of the two point correlation function and of the avalanche size distribution, one is confronted with the problem of estimating the minimal and maximal eigenvalues (denoted by λ − (M ) and λ + (M )) of a product of matrices of the form where the x i 's are the characteristic ratios of some recursively defined subtrees. When x i = 1 for all i, this exactly corresponds to the analysis in [5]   where o(1) tends to zero as n → ∞, uniformly in the choice of the x i 's.
PROOF. The eigenvalues are given by with a = Tr(M ) and b = det(M ) (we abbreviate M for M (x 1 , . . . , x n ) in these expressions). For n = 1, a 2 − 4b = 5 + 4x ≥ 0. Hence For n ≥ 2, we estimate, as in [11] a = Tr(M ) ≥ n i=1 (2 + x i ) and using In particular, 4b ≤ a 2 which implies that the eigenvalues are real and nonnegative (a ≥ 0). Inequality (31) then follows from (36), (33). Given that (by (36)) h = b/a 2 tends to zero as n → ∞ at a speed at least C(4/9) n , we have which proves the third statement.   where y i = (1 + x i ) and where The result then follows from expansion of this product. 1. For all n ≥ 1, 2. For all r ≥ 1, k 1 , . . . , k r , k r+1 ≥ 0, and PROOF. Identity (42) follows from (41) and invariance of the trace under cyclic permutations. To prove (41), use the expression in (40) for E n 1 E 2 , and estimate the diagonal elements of the product which implies the result. (1 + 2y i ) 2 −n(16/25) .
As a consequence we have the following uniform upper bound , -with the convention α n+1 = 1 -is the number of intervals of successive 1's in the configuration α. Hence by (38) we obtain the lower bound Tr (M (x 1 , . . . , x n ) defines a (y i ) i dependent probability measure on the α's. Now apply Jensen's inequality in (46) to obtain Finally, use Finally we use y i = (1 + x i ) with x i ∈ [1/2, 1] and the fact that x → (1 + x)(3 + 2x) −1 is increasing to estimate which implies (44).

Transfer matrix: annealed estimates
In this section we look at the eigenvalues of M (x 1 , . . . , x n ) where now the x i 's are i.i.d. with a law µ on [1/2, 1]. We denote by P the joint law of the x i 's and by E the corresponding expectation. We start with the following lemma.
Then we have: 1. Concentration property: there exists C > 0 such that for all ε > 0

The limits
exist and satisfy L ± = E(L ± ) almost surely; moreover
6 Covariance and avalanche sizes 6.1 Quenched and annealed covariance THEOREM 6.1. Let T be a finite or infinite subtree of the full binary tree B ∞ . As before, denote by µ T the uniform measure on the recurrent configurations R T of the sandpile model on T and let u, v ∈ T be at mutual distance n. Then we have the following estimate for the two point correlation function In particular, Cov(u, v, T ) is absolutely summable, uniformly in the choice of the subtree T : sup PROOF. The first statement (56) follows from the expression of the covariance in Subsection 2.5.3, eq. (5), together with Proposition 4.1. Since the number of vertices at distance n in T from the origin is bounded from above by 2 n , and γ < 1/2, (57) follows from (56). Then we have the following annealed lower bound on the covariance.
PROOF. The result follows from the expression of the covariance in Subsection 2.5.3, eq. (5), together with (52).

Avalanche sizes
We start by defining a matrix associated to an avalanche cluster. Roughly speaking, the probability that the avalanche coincides with the cluster is the inverse of the maximal eigenvalue of this matrix. Let T ⊂ B ∞ be a finite or infinite subtree of the full rootless binary tree, containing the origin o. Let C be a connected subset of T containing the origin. DEFINITION 6.1. Let |C | = n. The matrix associated to C in T , denoted by M (C ) is defined as follows. To the vertices in C are associated n + 2 subtrees T 1 , . . . , T n+2 with corresponding characteristic ratios x(T i ). Then (59) LEMMA 6.1. As before, for a stable height configuration η, let Av(o, η) denote the avalanche caused by addition of a single grain at the origin in η, and µ T the uniform measure on recurrent configurations R T on T . Then there exist constants c 1 , c 2 > 0 such that where M (C ) is the matrix of (59) associated to C of Definition 6.1.
PROOF. It follows from the expression for µ T (Av(o, η) = C ) in Subsection 2.5.3, eq. (6). The growth rate gives the dominant exponential factor in the number A n (o, T ). E.g. if T is the full binary tree, κ = log 4, since see also [5], Subsection 2.5.2. For a stationary binary branching process we have the upper bound For an exact expression E(A n (o, T )) from which inequality (63) follows immediately, we refer to the Appendix, Proposition 7.2.
In particular if the tree T has growth rate κ(T ) < 34 25 log 2, then the avalanche size decays exponentially.
2. For the stationary binary branching tree with branching probability p, we have the estimate (65) In particular, for p < 0.54511...
the averaged (over the realization of the tree) probability of an avalanche size n decays exponentially in n.
the averaged (over the realization of the tree) probability of an avalanche size n decays exponentially in n.
PROOF. The first result, i.e., (64), follows from Lemma 6.1, and the estimate (43) which gives the following lower bound on the largest eigenvalue The second result follows from Lemma 6.1, (67), (63), (62) and Proposition 7.1 below. The solution of is given by The l.h.s. of (68) as a function of p is monotone, which yields that for p < 0.54511 the expected number of avalanches of size n decays exponentially. The third statement follows from Lemma 6.1, (67), (63), (62) and Remark 7.1. REMARK 6.1. From Theorem 6.3 we conclude that avalanche sizes decay exponentially for p small enough. We know that on the full binary tree (corresponding to p = 1) avalanche sizes have a power law decay. It is a natural question whether the transition between exponential and power law decay occurs at some unique non-trivial value p c ∈ (1/2, 1). This question is related to the behavior of a a random walk on the full binary tree, killed upon exiting a random subtree. Would this random walk have a survival time with a finite exponential moment, then we would have p c = 1. Notice that for the corresponding problem on the lattice Z d , this random walk does not have a finite exponential moment of its survival time, because the tail of the survival time is stretched exponential (Donsker-Varadhan tail e −ct d/(d+2) ). However, on the tree in the annealed case corresponding to Theorem 6.3, the decay of the tail of the survival time is not straightforward. Formally taking the limit d → ∞ in the Donsker-Varadhan tail suggests that on the tree, the survival time of this random walk has a finite exponential moment, which points into the direction p c = 1. See also [2] for trapped random walk on a tree. which follows from the observation that the number of n vertex animals containing the root of the full binary tree is bounded from above by 4 n , and in the case of binary branching tree, each of these vertices is present with probability p independently, which gives the factor p n .
PROOF. If we define a n to be the expected number of connected clusters containing o and with n vertices, then a n satisfies the recursion relation a n = 1l {n=1} + p1l {n≥2} n−1 i=0 a i a n−i−1 , n ≥ 1.
Furthermore, by definition we put a 0 = 1. Thus, going from vertices to edges, Introduce then the associated generating function: A(x) = ∞ n=0 a n x n = 1 + x + p ∞ n=2 n−1 i=0 a i a n−i−1 x n = 1 + x + px(A 2 (x) − 1). For , this power series is convergent and, using A(0) = 1, it is given by . Use (for z such that 4z < 1) to obtain We put k = n + j, then n = k − j and j ≤ n + 1 yields j ≤ (k + 1)/2 , hence the expected number of clusters a k containing the origin is equal to where b j,k = 2(k − j) k − j k − j + 1 j 1 − p p with c = (1 − p)/p. All terms in the sum (78) are exponentially large in k, hence the exponential growth of the sum is determined by its maximal term. Using Stirling's approximation n! ≈ n n e −n √ 2πn for the right hand side of (79), we obtain the upper bound Define, for x ∈ [0, (k + 1)/2], This function ϕ attains its maximum at x = k(1 − √ p)/2, which combined with (80), implies Plugging (81) into (78) and bounding from above induces (70).
In the following proposition we give an exact formula for E(A n (o, T o )).