Strong Spherical Asymptotics for Rotor-Router
Aggregation and the Divisible Sandpile
Lionel Levine∗ and Yuval Peres†
University of California, Berkeley and Microsoft Research
October 7, 2008
Abstract
The rotor-router model is a deterministic analogue of random walk.
It can be used to define a deterministic growth model analogous to
internal DLA. We prove that the asymptotic shape of this model is a
Euclidean ball, in a sense which is stronger than our earlier work [14].
For the shape consisting of n = ωdrd sites, where ωd is the volume of
the unit ball in Rd, we show that the inradius of the set of occupied
sites is at least r −O(log r), while the outradius is at most r +O(rα)
for any α > 1 − 1/d. For a related model, the divisible sandpile, we
show that the domain of occupied sites is a Euclidean ball with error in
the radius a constant independent of the total mass. For the classical
abelian sandpile model in two dimensions, with n = πr2 particles, we
show that the inradius is at least r/
3, and the outradius is at most
(r + o(r))/
2. This improves on bounds of Le Borgne and Rossin.
Similar bounds apply in higher dimensions, improving on bounds of
Fey and Redig.
1 Introduction
Rotor-router walk is a deterministic analogue of random walk, first studied
by Priezzhev et al. [18] under the name “Eulerian walkers.” At each site
∗supported by an NSF Graduate Research Fellowship, and NSF grant DMS-0605166;
levine(at)math.berkeley.edu
†partially supported by NSF grant DMS-0605166; peres(at)stat.berkeley.edu
Keywords: abelian sandpile, asymptotic shape, discrete Laplacian, divisible sand-
pile, growth model, internal diffusion limited aggregation, rotor-router model
2000 Mathematics Subject Classifications: Primary 60G50; Secondary 35R35
Figure 1: Rotor-router aggregate of one million particles in Z2. Each site is
colored according to the direction of its rotor.
in the integer lattice Z2 is a rotor pointing north, south, east or west. A
particle starts at the origin; during each time step, the rotor at the particle’s
current location is rotated clockwise by 90 degrees, and the particle takes a
step in the direction of the newly rotated rotor. In rotor-router aggregation,
introduced by Jim Propp, we start with n particles at the origin; each par-
ticle in turn performs rotor-router walk until it reaches a site not occupied
by any other particles. Let An denote the resulting region of n occupied
sites. For example, if all rotors initially point north, the sequence will begin
A1 = {(0, 0)}, A2 = {(0, 0), (1, 0)}, A3 = {(0, 0), (1, 0), (0,−1)}. The region
A1,000,000 is pictured in Figure 1. In higher dimensions, the model can be
defined analogously by repeatedly cycling the rotors through an ordering of
the 2d cardinal directions in Zd.
Jim Propp observed from simulations in two dimensions that the regions
An are extraordinarily close to circular, and asked why this was so [7, 19].
Despite the impressive empirical evidence for circularity, the best result
known until now [14] says only that if An is rescaled to have unit volume,
the volume of the symmetric difference of An with a ball of unit volume tends
to zero as a power of n, as n ↑ ∞. The main outline of the argument is
summarized in [15]. Fey and Redig [5] also show that An contains a diamond.
In particular, these results do not rule out the possibility of “holes” in An
far from the boundary or of long tendrils extending far beyond the boundary
of the ball, provided the volume of these features is negligible compared to
Our main result is the following, which rules out the possibility of holes
far from the boundary or of long tendrils in the rotor-router shape. For
r ≥ 0 let
Br = {x ∈ Zd : |x| < r}.
Theorem 1.1. Let An be the region formed by rotor-router aggregation in
Zd starting from n particles at the origin and any initial rotor state. There
exist constants c, c′ depending only on d, such that
Br−c log r ⊂ An ⊂ Br(1+c′r−1/d log r)
where r = (n/ωd)1/d, and ωd is the volume of the unit ball in Rd.
We remark that the same result holds when the rotors evolve according
to stacks of bounded discrepancy; see the remark following Lemma 5.1.
Internal diffusion limited aggregation (“internal DLA”) is an analogous
growth model defined using random walks instead of rotor-router walks.
Starting with n particles at the origin, each particle in turn performs simple
random walk until it reaches an unoccupied site. Lawler, Bramson and
Griffeath [10] showed that for internal DLA in Zd, the occupied region An,
rescaled by a factor of n1/d, converges with probability one to a Euclidean
ball in Rd as n→∞. Lawler [11] estimated the rate of convergence. By way
of comparison with Theorem 1.1, if In is the internal DLA region formed
from n particles started at the origin, the best known bounds [11] are (up
to logarithmic factors)
Br−r1/3 ⊂ In ⊂ Br+r1/3
for all sufficiently large n, with probability one.
We also study another model which is slightly more difficult to define,
but much easier to analyze. In the divisible sandpile, each site x ∈ Zd
starts with a quantity of “mass” ν0(x) ∈ R≥0. A site topples by keeping up
Figure 2: Classical abelian sandpile aggregate of one million particles in Z2.
Colors represent the number of grains at each site.
to mass 1 for itself, and distributing the excess (if any) equally among its
neighbors. Thus if x has mass m > 1, then each of the 2d neighboring sites
gains mass (m−1)/2d when we topple x, and x is left with mass 1; if m ≤ 1,
then no mass moves when we topple x.
Note that individual topplings do not commute; however, the divisible
sandpile is “abelian” in the following sense.
Proposition 1.2. Let x1, x2, . . . ∈ Zd be a sequence with the property that
for any x ∈ Zd there are infinitely many terms xk = x. Let
uk(x) = total mass emitted by x after toppling x1, . . . , xk;
νk(x) = amount of mass present at x after toppling x1, . . . , xk.
Then uk ↑ u and νk → ν ≤ 1. Moreover, the limits u and ν are independent
of the sequence {xk}.
The abelian property can be generalized as follows: after performing
some topplings, we can add some additional mass and then continue top-
pling. The resulting limits u and ν will be the same as in the case when all
mass was initially present. For a further generalization, see [16].
The limiting function u in Proposition 1.2 is the odometer function for
the divisible sandpile. This function plays a central role in our analysis. The
limit ν represents the final mass distribution. Sites x ∈ Zd with ν(x) = 1
are called fully occupied. Proposition 1.2 is proved in section 3, along with
the following.
Theorem 1.3. For m ≥ 0 let Dm ⊂ Zd be the domain of fully occupied
sites for the divisible sandpile formed from a pile of mass m at the origin.
There exist constants c, c′ depending only on d, such that
Br−c ⊂ Dm ⊂ Br+c′ ,
where r = (m/ωd)1/d and ωd is the volume of the unit ball in Rd.
The divisible sandpile is similar to the “oil game” studied by Van den
Heuvel [22]. In the terminology of [5], it also corresponds to the h → −∞
limit of the classical abelian sandpile (defined below), that is, the abelian
sandpile started from the initial condition in which every site has a very
deep “hole.”
In the classical abelian sandpile model [1], each site in Zd has an integer
number of grains of sand; if a site has at least 2d grains, it topples, sending
one grain to each neighbor. If n grains of sand are started at the origin in
Zd, write Sn for the set of sites that are visited during the toppling process;
in particular, although a site may be empty in the final state, we include it
in Sn if it was occupied at any time during the evolution to the final state.
Until now the best known constraints on the shape of Sn in two dimen-
sions were due to Le Borgne and Rossin [12], who proved that
{x ∈ Z2 |x1 + x2 ≤
n/12− 1} ⊂ Sn ⊂ {x ∈ Z2 |x1, x2 ≤
n/2}.
Fey and Redig [5] proved analogous bounds in higher dimensions, and ex-
tended these bounds to arbitrary values of the height parameter h. This
parameter is discussed in section 4.
The methods used to prove the near-perfect circularity of the divisible
sandpile shape in Theorem 1.3 can be used to give constraints on the shape
of the classical abelian sandpile, improving on the bounds of [5] and [12].
Theorem 1.4. Let Sn be the set of sites that are visited by the classical
abelian sandpile model in Zd, starting from n particles at the origin. Write
n = ωdrd. Then for any � > 0 we have
Bc1r−c2 ⊂ Sn ⊂ Bc′1r+c′2
Figure 3: Known bounds on the shape of the classical abelian sandpile in
Z2. The inner diamond and outer square are due to Le Borgne and Rossin
[12]; the inner and outer circles are those in Theorem 1.4.
where
c1 = (2d− 1)−1/d, c′1 = (d− �)
−1/d.
The constant c2 depends only on d, while c′2 depends only on d and �.
Note that Theorem 1.4 does not settle the question of the asymptotic
shape of Sn, and indeed it is not clear from simulations whether the asymp-
totic shape in two dimensions is a disc or perhaps a polygon (Figure 2). To
our knowledge even the existence of an asymptotic shape is not known.
The rest of the paper is organized as follows. In section 2, we derive the
basic Green’s function estimates that are used in the proofs of Theorems 1.1,
1.3 and 1.4. In section 3 we prove Proposition 1.2 and Theorem 1.3 for the
divisible sandpile. In section 4 we adapt the methods of the previous section
to prove Theorem 1.4 for the classical abelian sandpile model. Section 5 is
devoted to the proof of Theorem 1.1.
2 Basic Estimate
Write (Xk)k≥0 for simple random walk in Zd, and for d ≥ 3 denote by
g(x) = Eo#{k|Xk = x}
the expected number of visits to x by simple random walk started at the
origin. This is the discrete harmonic Green’s function in Zd; it satisfies
∆g(x) = 0 for x 6= o, and ∆g(o) = −1, where ∆ is the discrete Laplacian
∆g(x) =
g(y)− g(x).
The sum is over the 2d lattice neighbors y of x. In dimension d = 2, simple
random walk is recurrent, so the expectation defining g(x) is infinite. Here
we define the potential kernel
g(x) = lim
gn(x)− gn(o) (1)
where
gn(x) = Eo#{k ≤ n|Xk = x}.
The limit defining g(x) in (1) is finite [9, 20], and g(x) has Laplacian ∆g(x) =
0 for x 6= o, and ∆g(o) = −1. Note that (1) is the negative of the usual
definition of the potential kernel; we have chosen this sign convention so
that g has the same Laplacian in dimension two as in higher dimensions.
Fix a real number m > 0 and consider the function on Zd
γ̃d(x) = |x|2 +mg(x). (2)
Let r be such that m = ωdrd, and let
γd(x) = γ̃d(x)− γ̃d(brce1) (3)
where e1 is the first standard basis vector in Zd. The function γd plays a
central role in our analysis. To see where it comes from, recall the divisible
sandpile odometer function of Proposition 1.2
u(x) = total mass emitted from x.
Let Dm ⊂ Zd be the domain of fully occupied sites for the divisible sandpile
formed from a pile of mass m at the origin. For x ∈ Dm, since each neighbor
y of x emits an equal amount of mass to each of its 2d neighbors, we have
∆u(x) =
u(y)− u(x)
= mass received by x−mass emitted by x
= 1−mδox.
Moreover, u = 0 on ∂Dm. By construction, the function γd obeys the same
Laplacian condition: ∆γd = 1 −mδo; and as we will see shortly, γd ≈ 0 on
∂Br. Since we expect the domain Dm to be close to the ball Br, we should
expect that u ≈ γd. In fact, we will first show that u is close to γd, and then
use this to conclude that Dm is close to Br.
We will use the following estimates for the Green’s function [6, 21]; see
also [9, Theorems 1.5.4 and 1.6.2].
g(x) =
log |x|+ κ+O(|x|−2), d = 2
ad|x|2−d +O(|x|−d), d ≥ 3.
Here ad =
(d−2)ωd
, where ωd is the volume of the unit ball in Rd, and κ is
a constant whose value we will not need to know. For x ∈ Zd we use |x| to
denote the Euclidean norm of x. Here and throughout the paper, constants
in error terms denoted O(·) depend only on d.
We will need an estimate for γd near the boundary of the ball Br. We
first consider dimension d = 2. From (4) we have
γ̃2(x) = φ(x)− κm+O(m|x|−2), (5)
where
φ(x) = |x|2 −
log |x|.
In the Taylor expansion of φ about |x| = r
φ(x) = φ(r)− φ′(r)(r − |x|) +
φ′′(t)(r − |x|)2 (6)
the linear term vanishes, leaving
γ2(x) =
(r − |x|)2 +O(m|x|−2) (7)
for some t between |x| and r.
In dimensions d ≥ 3, from (4) we have
γ̃d(x) = |x|2 + adm|x|2−d +O(m|x|−d).
Setting φ(x) = |x|2 +adm|x|2−d, the linear term in the Taylor expansion (6)
of φ about |x| = r again vanishes, yielding
γd(x) =
1 + (d− 1)(r/t)d
(r − |x|)2 +O(m|x|−d)
for t between |x| and r. Together with (7), this yields the following estimates
in all dimensions d ≥ 2.
Lemma 2.1. Let γd be given by (3). For all x ∈ Zd we have
γd(x) ≥ (r − |x|)2 +O
. (8)
Lemma 2.2. Let γd be given by (3). Then uniformly in r,
γd(x) = O(1), x ∈ Br+1 −Br−1.
The following lemma is useful for x near the origin, where the error term
in (8) blows up.
Lemma 2.3. Let γd be given by (3). Then for sufficiently large r, we have
γd(x) >
, ∀x ∈ Br/3.
Proof. Since γd(x) − |x|2 is superharmonic, it attains its minimum in Br/3
at a point z on the boundary. Thus for any x ∈ Br/3
γd(x)− |x|2 ≥ γd(z)− |z|2,
hence by Lemma 2.1
γd(x) ≥ (2r/3)2 − (r/3)2 +O(1) >
Lemmas 2.1 and 2.3 together imply the following.
Lemma 2.4. Let γd be given by (3). There is a constant a depending only
on d, such that γd ≥ −a everywhere.
3 Divisible Sandpile
Let µ0 be a nonnegative function on Zd with finite support. We start with
mass µ0(y) at each site y. The operation of toppling a vertex x yields the
mass distribution
Txµ0 = µ0 + α(x)∆δx
where α(x) = max(µ0(x)−1, 0) and ∆ is the discrete Laplacian on Zd. Thus
if µ0(x) ≤ 1 then Txµ0 = µ0 and no mass topples; if µ0(x) > 1 then the
mass in excess of 1 is distributed equally among the neighbors of x.
Let x1, x2, . . . ∈ Zd be a sequence with the property that for any x ∈ Zd
there are infinitely many terms xk = x. Let
µk(y) = Txk . . . Tx1µ0(y).
be the amount of mass present at y after toppling the sites x1, . . . , xk in
succession. The total mass emitted from y during this process is
uk(y) :=
j≤k:xj=y
µj−1(y)− µj(y) =
j≤k:xj=y
αj(y) (9)
where αj(y) = max(µj(y)− 1, 0).
Lemma 3.1. As k ↑ ∞ the functions uk and µk tend to limits uk ↑ u and
µk → µ. Moreover, these limits satisfy
µ = µ0 + ∆u ≤ 1.
Proof. Write M =
x µ0(x) for the total starting mass, and let B ⊂ Z
be a ball centered at the origin containing all points within L1-distance M
of the support of µ0. Note that if µk(x) > 0 and µ0(x) = 0, then x must
have received its mass from a neighbor, so µk(y) ≥ 1 for some y ∼ x. Since∑
x µk(x) = M , it follows that µk is supported on B. Let R be the radius
of B, and consider the quadratic weight
µk(x)|x|2 ≤MR2.
Since µk(xk)−µk−1(xk) = −αk(xk) and for y ∼ xk we have µk(y)−µk−1(y) =
αk(xk), we obtain
Qk −Qk−1 = αk(xk)
|y|2 − |xk|2
= αk(xk).
Summing over k we obtain from (9)
Qk = Q0 +
uk(x).
Fixing x, the sequence uk(x) is thus increasing and bounded above, hence
convergent.
Given neighboring vertices x ∼ y, since y emits an equal amount of mass
to each of its 2d neighbors, it emits mass uk(y)/2d to x up to time k. Thus
x receives a total mass of 1
y∼x uk(y) from its neighbors up to time k.
Comparing the amount of mass present at x before and after toppling, we
obtain
µk(x) = µ0(x) + ∆uk(x).
Since uk ↑ u we infer that µk → µ := µ0 + ∆u. Note that if xk = x, then
µk(x) ≤ 1. Since for each x ∈ Zd this holds for infinitely many values of k,
the limit satisfies µ ≤ 1.
A function s on Zd is superharmonic if ∆s ≤ 0. Given a function γ on
Zd the least superharmonic majorant of γ is the function
s(x) = inf{f(x) | f is superharmonic and f ≥ γ}.
The study of the least superharmonic majorant is a classical topic in analysis
and PDE; see, for example, [8]. Note that if f is superharmonic and f ≥ γ
f(x) ≥
f(y) ≥
s(y).
Taking the infimum on the left side we obtain that s is superharmonic.
Lemma 3.2. The limit u in Lemma 3.1 is given by u = s+ γ, where
γ(x) = |x|2 +
g(x− y)µ0(y)
and s is the least superharmonic majorant of −γ.
Proof. By Lemma 3.1 we have
∆u = µ− µ0 ≤ 1− µ0.
Since ∆γ = 1−µ0, the difference u−γ is superharmonic. As u is nonnegative,
it follows that u − γ ≥ s. For the reverse inequality, note that s + γ − u is
superharmonic on the domain D = {x | µ(x) = 1} of fully occupied sites
and is nonnegative outside D, hence nonnegative inside D as well.
As a corollary of Lemmas 3.1 and 3.2, we obtain the abelian property of
the divisible sandpile, Proposition 1.2, which was stated in the introduction.
We now turn to the case of a point source mass m started at the origin:
µ0 = mδo. More general starting distributions are treated in [16], where
we identify the scaling limit of the divisible sandpile model and show that
it coincides with that of internal DLA and of the rotor-router model. In
the case of a point source of mass m, the natural question is to identify
the shape of the resulting domain Dm of fully occupied sites, i.e. sites x for
which µ(x) = 1. According to Theorem 3.3, Dm is extremely close to a ball
of volume m; in fact, the error in the radius is a constant independent of m.
As before, for r ≥ 0 we write
Br = {x ∈ Zd : |x| < r}
for the lattice ball of radius r centered at the origin.
Theorem 3.3. For m ≥ 0 let Dm ⊂ Zd be the domain of fully occupied sites
for the divisible sandpile formed from a pile of size m at the origin. There
exist constants c, c′ depending only on d, such that
Br−c ⊂ Dm ⊂ Br+c′ ,
where r = (m/ωd)1/d and ωd is the volume of the unit ball in Rd.
The idea of the proof is to use Lemma 3.2 along with the basic estimates
on γ, Lemmas 2.1 and 2.2, to obtain estimates on the odometer function
u(x) = total mass emitted from x.
We will need the following simple observation.
Lemma 3.4. For every point x ∈ Dm − {o} there is a path x = x0 ∼ x1 ∼
. . . ∼ xk = o in Dm with u(xi+1) ≥ u(xi) + 1.
Proof. If xi ∈ Dm − {o}, let xi+1 be a neighbor of xi maximizing u(xi+1).
Then xi+1 ∈ Dm and
u(xi+1) ≥
= u(xi) + ∆u(xi)
= u(xi) + 1,
where in the last step we have used the fact that xi ∈ Dm.
Proof of Theorem 3.3. We first treat the inner estimate. Let γd be given by
(3). By Lemma 3.2 the function u−γd is superharmonic, so its minimum in
the ball Br is attained on the boundary. Since u ≥ 0, we have by Lemma 2.2
u(x)− γd(x) ≥ −C, x ∈ ∂Br
for a constant C depending only on d. Hence by Lemma 2.1,
u(x) ≥ (r − |x|)2 − C ′rd/|x|d, x ∈ Br. (10)
for a constant C ′ depending only on d. It follows that there is a constant c,
depending only on d, such that u(x) > 0 whenever r/3 ≤ |x| < r − c. Thus
Br−c−Br/3 ⊂ Dm. For x ∈ Br/3, by Lemma 2.3 we have u(x) ≥ r2/4−C >
0, hence Br/3 ⊂ Dm.
For the outer estimate, note that u−γd is harmonic onDm. By Lemma 2.4
we have γd ≥ −a everywhere, where a depends only on d. Since u vanishes
on ∂Dm it follows that u − γd ≤ a on Dm. Now for any x ∈ Dm with
r − 1 < |x| ≤ r, we have by Lemma 2.2
u(x) ≤ γd(x) + a ≤ c′
for a constant c′ depending only on d. Lemma 3.4 now implies that Dm ⊂
Br+c′+1.
4 Classical Sandpile
We consider a generalization of the classical abelian sandpile, proposed by
Fey and Redig [5]. Each site in Zd begins with a “hole” of depth H. Thus,
each site absorbs the first H grains it receives, and thereafter functions
normally, toppling once for each additional 2d grains it receives. If H is
negative, we can interpret this as saying that every site starts with h = −H
grains of sand already present. Aggregation is only well-defined in the regime
h ≤ 2d− 2, since for h = 2d− 1 the addition of a single grain already causes
every site in Zd to topple infinitely often.
Let Sn,H be the set of sites that are visited if n particles start at the
origin in Zd. Fey and Redig [5, Theorem 4.7] prove that
lim sup
# (Sn,H 4BH−1/dr) = 0,
where n = ωdrd, and 4 denotes symmetric difference. The following theo-
rem strengthens this result.
Theorem 4.1. Fix an integer H ≥ 2−2d. Let Sn = Sn,H be the set of sites
that are visited by the classical abelian sandpile model in Zd, starting from
n particles at the origin, if every lattice site begins with a hole of depth H.
Write n = ωdrd. Then
Bc1r−c2 ⊂ Sn,H
where
c1 = (2d− 1 +H)−1/d
and c2 is a constant depending only on d. Moreover if H ≥ 1− d, then for
any � > 0 we have
Sn,H ⊂ Bc′1r+c′2
where
c′1 = (d− �+H)
and c′2 is independent of n but may depend on d, H and �.
Note that the ratio c1/c′1 ↑ 1 as H ↑ ∞. Thus, the classical abelian
sandpile run from an initial state in which each lattice site starts with a
deep hole yields a shape very close to a ball. Intuitively, one can think of the
classical sandpile with deep holes as approximating the divisible sandpile,
whose limiting shape is a ball by Theorem 3.3. Following this intuition,
we can adapt the proof of Theorem 3.3 to prove Theorem 4.1; just one
additional averaging trick is needed, which we explain below.
Consider the odometer function for the abelian sandpile
u(x) = total number of grains emitted from x.
Let Tn = {x|u(x) > 0} be the set of sites which topple at least once. Then
Tn ⊂ Sn ⊂ Tn ∪ ∂Tn.
In the final state, each site which has toppled retains between 0 and 2d− 1
grains, in addition to the H that it absorbed. Hence
H ≤ ∆u(x) + nδox ≤ 2d− 1 +H, x ∈ Tn. (11)
We can improve the lower bound by averaging over a small box. For x ∈ Zd
Qk(x) = {y ∈ Zd : ||x− y||∞ ≤ k}
be the box of side length 2k + 1 centered at x, and let
u(k)(x) = (2k + 1)−d
y∈Qk(x)
u(y).
Write
T (k)n = {x |Qk(x) ⊂ Tn}.
Le Borgne and Rossin [12] observe that if T is a set of sites all of which
topple, the number of grains remaining in T is at least the number of edges
internal to T : indeed, for each internal edge, the endpoint that topples last
sends the other a grain which never moves again. Since the box Qk(x) has
2dk(2k + 1)d−1 internal edges, we have
∆u(k)(x) ≥
2k + 1
d+H −
(2k + 1)d
1Qk(o)(x), x ∈ T
n . (12)
The following lemma is analogous to Lemma 3.4.
Lemma 4.2. For every point x ∈ Tn adjacent to ∂Tn there is a path x =
x0 ∼ x1 ∼ . . . ∼ xm = o in Tn with u(xi+1) ≥ u(xi) + 1.
Proof. By (11) we have
u(y) ≥ u(xi).
Since u(xi−1) < u(xi), some term u(y) in the sum above must exceed u(xi).
Let xi+1 = y.
Proof of Theorem 4.1. Let
ξ̃d(x) = (2d− 1 +H)|x|2 + ng(x),
and let
ξd(x) = ξ̃d(x)− ξ̃d(bc1rce1).
Taking m = n/(2d− 1 +H) in Lemma 2.2, we have
u(x)− ξd(x) ≥ −ξd(x) ≥ −C(2d− 1 +H), x ∈ ∂Bc1r (13)
for a constant C depending only on d. By (11), u− ξd is superharmonic, so
u − ξd ≥ −C(2d − 1 + H) in all of Bc1r. Hence by Lemma 2.1 we have for
x ∈ Bc1r
u(x) ≥ (2d− 1 +H)
(c1r − |x|)2 − C ′(c1r)d/|x|d
, (14)
where C ′ depends only on d. It follows that u is positive on Bc1r−c2−Bc1r/3
for a suitable constant c2 depending only on d. For x ∈ Bc1r/3, by Lemma 2.3
we have u(x) > (2d− 1 +H)(c21r
2/4− C) > 0. Thus Bc1r−c2 ⊂ Tn ⊂ Sn.
For the outer estimate, let
ψ̂d(x) = (d− �+H)|x|2 + ng(x).
Choose k large enough so that 2k
d ≥ d− �, and define
ψ̃d(x) = (2k + 1)
y∈Qk(x)
ψ̂d(y).
Finally, let
ψd(x) = ψ̃d(x)− ψ̃d(bc′1rce1).
By (12), u(k) − ψd is subharmonic on T
n . Taking m = n/(d − � + H)
in Lemma 2.4, there is a constant a depending only on d, such that ψd ≥
−a(d+H) everywhere. Since u(k) ≤ (2d+H)(d+1)k on ∂T (k)n it follows that
u(k) − ψd ≤ a(d + H) + (2d + H)(d+1)k on T
n . Now for any x ∈ Sn with
c′1r − 1 < |x| ≤ c
1r we have by Lemma 2.2
u(k)(x) ≤ ψd(x) + a(d+H) + (2d+H)(d+1)k ≤ c̃2
for a constant c̃2 depending only on d, H and �. Then u(x) ≤ c′2 := (2k +
1)dc̃2. Lemma 4.2 now implies that Tn ⊂ Bc′1r+c′2 , and hence
Sn ⊂ Tn ∪ ∂Tn ⊂ Bc′1r+c′2+1.
We remark that the crude bound of (2d + H)(d+1)k used in the proof
of the outer estimate can be improved to a bound of order k2H, and the
final factor of (2k + 1)d can be replaced by a constant factor independent
of k and H, using the fact that a nonnegative function on Zd with bounded
Laplacian cannot grow faster than quadratically; see [16].
5 Rotor-Router Model
Given a function f on Zd, for a directed edge (x, y) write
∇f(x, y) = f(y)− f(x).
Given a function s on directed edges in Zd, write
div s(x) =
s(x, y).
The discrete Laplacian of f is then given by
∆f(x) = div∇f =
f(y)− f(x).
5.1 Inner Estimate
Fixing n ≥ 1, consider the odometer function for rotor-router aggregation
u(x) = total number exits from x by the first n particles.
We learned the idea of using the odometer function to study the rotor-router
shape from Matt Cook [2].
Lemma 5.1. For a directed edge (x, y) in Zd, denote by κ(x, y) the net
number of crossings from x to y performed by the first n particles in rotor-
router aggregation. Then
∇u(x, y) = −2dκ(x, y) +R(x, y) (15)
for some edge function R which satisfies
|R(x, y)| ≤ 4d− 2
for all edges (x, y).
Remark. In the more general setting of rotor stacks of bounded discrepancy,
the 4d− 2 will be replaced by a different constant here.
Proof. Writing N(x, y) for the number of particles routed from x to y, we
u(x)− 2d+ 1
≤ N(x, y) ≤
u(x) + 2d− 1
hence
|∇u(x, y) + 2dκ(x, y)| = |u(y)− u(x) + 2dN(x, y)− 2dN(y, x)|
≤ 4d− 2.
In what follows, C0, C1, . . . denote constants depending only on d.
Lemma 5.2. Let Ω ⊂ Zd − {o} with 2 ≤ #Ω <∞. Then∑
|y|1−d ≤ C0 Diam(Ω).
Proof. For each positive integer k, let
Sk = {y ∈ Zd : k ≤ |y| < k + 1}.
Then ∑
|y|1−d ≤ k1−d#Sk ≤ C ′0
for a constant C ′0 depending only on d. Since Ω can intersect at most
Diam(Ω) + 1 ≤ 2 Diam(Ω) distinct sets Sk, taking C0 = 2C ′0 the proof is
complete.
Lemma 5.3. Let G = GBr be the Green’s function for simple random walk
in Zd stopped on exiting Br. For any ρ ≥ 1 and x ∈ Br,∑
|x−y|≤ρ
|G(x, y)−G(x, z)| ≤ C1ρ. (16)
Proof. Let (Xt)t≥0 denote simple random walk in Zd, and let T be the first
exit time from Br. For fixed y, the function
A(x) = g(x− y)− Exg(XT − y) (17)
has Laplacian ∆A(x) = −δxy in Br and vanishes on ∂Br, hence A(x) =
G(x, y).
Let x, y ∈ Br and z ∼ y. From (4) we have
|g(x− y)− g(x− z)| ≤
|x− y|d−1
, y, z 6= x.
Using the triangle inequality together with (17), we obtain
|G(x, y)−G(x, z)| ≤ |g(x− y)− g(x− z)|+ Ex|g(XT − y)− g(XT − z)|
|x− y|d−1
w∈∂Br
Hx(w)
|w − y|d−1
where Hx(w) = Px(XT = w).
Write D = {y ∈ Br : |x− y| ≤}. Then∑
y 6=x
z 6=x
|G(x, y)−G(x, z)| ≤ C3ρ+ C2
w∈∂Br
Hx(w)
|w − y|1−d. (18)
Figure 4: Diagram for the Proof of Lemma 5.4.
Taking Ω = w − D in Lemma 5.2, the inner sum on the right is at most
C0Diam(D) ≤ 2C0ρ, so the right side of (18) is bounded above by C1ρ for
a suitable C1.
Finally, the terms in which y or z coincides with x make a negligible
contribution to the sum in (16), since for y ∼ x ∈ Zd
|G(x, x)−G(x, y)| ≤ |g(o)−g(x−y)|+Ex|g(XT −x)−g(XT −y)| ≤ C4.
Lemma 5.4. Let H1, H2 be linear half-spaces in Zd, not necessarily parallel
to the coordinate axes. Let Ti be the first hitting time of Hi. If x /∈ H1∪H2,
Px(T1 > T2) ≤
h1 + 1
where hi is the distance from x to Hi.
Proof. If one of H1, H2 contains the other, the result is vacuous. Otherwise,
let H̃i be the half-space shifted parallel to Hci by distance 2h2 in the direction
of x, and let T̃i be the first hitting time of Hi∪H̃i. Let (Xt)t≥0 denote simple
random walk in Zd, and write Mt for the (signed) distance from Xt to the
hyperplane defining the boundary of H1, with M0 = h1. Then Mt is a
martingale with bounded increments. Since ExT̃1 < ∞, we obtain from
optional stopping
h1 = ExMeT1 ≥ 2h2 Px (XeT1 ∈ H̃1)− Px(XeT1 ∈ H1),
hence
Px (XeT1 ∈ H̃1) ≤ h1 + 12h2 . (19)
Likewise, dM2t − t is a martingale with bounded increments, giving
ExT̃1 ≤ d ExM2eT1
≤ d(2h2 + 1)2 Px (XeT1 ∈ H̃1)
≤ d(h1 + 1)(2h2 + 1)
. (20)
Let T = min(T̃1, T̃2). Denoting by Dt the distance from Xt to the
hyperplane defining the boundary of H2, the quantity
D2t + (2h2 −Dt)
is a martingale. Writing p = Px(T = T̃2) we have
dh22 = EN0 = ENT ≥ p
(2h2)
2 + (1− p)dh22 − ExT
≥ (1 + p)dh22 − ExT
hence by (20)
h1 + 1
Finally by (19)
P(T1 > T2) ≤ p+ P(XeT1 ∈ H̃1) ≤ 52
h1 + 1
Lemma 5.5. Let x ∈ Br and let ρ = r + 1− |x|. Let
S∗k = {y ∈ Br : 2
kρ < |x− y| ≤ 2k+1ρ}. (21)
Let τk be the first hitting time of S∗k , and T the first exit time from Br. Then
Px(τk < T ) ≤ C22−k.
Figure 5: Diagram for the proof of Lemma 5.5.
Proof. Let H be the outer half-space tangent to Br at the point z ∈ ∂Br
closest to x. Let Q be the cube of side length 2kρ/
d centered at x. Then
Q is disjoint from S∗k , hence
Px(τk < T ) ≤ Px(T∂Q < T ) ≤ Px(T∂Q < TH)
where T∂Q and TH are the first hitting times of ∂Q and H. Let H1, . . . ,H2d
be the outer half-spaces defining the faces of Q, so that Q = Hc1 ∩ . . .∩H
By Lemma 5.4 we have
Px(T∂Q < TH) ≤
Px(THi < TH)
dist(x,H) + 1
dist(x,Hi)
2 dist(x,Hi)
Since dist(x,H) = |x − z| ≤ ρ and dist(x,Hi) = 2k−1ρ/
d, and ρ ≥ 1,
taking C2 = 20 d3/2(1 +
d)2 completes the proof.
Lemma 5.6. Let G = GBr be the Green’s function for random walk stopped
on exiting Br. Let x ∈ Br and let ρ = r + 1− |x|. Then∑
|G(x, y)−G(x, z)| ≤ C3ρ log
Proof. Let S∗k be given by (21), and let
W = {w ∈ ∂(S∗k ∪ ∂S
k) : |w − x| < 2
be the portion of the boundary of the enlarged spherical shell S∗k ∪∂S
k lying
closer to x. Let τW be the first hitting time of W , and T the first exit time
from Br. For w ∈W let
Hx(w) = Px(XτW∧T = w).
For any y ∈ S∗k and z ∼ y, simple random walk started at x must hit W
before hitting either y or z, hence
|G(x, y)−G(x, z)| ≤
Hx(w)|G(w, y)−G(w, z)|.
For any y ∈ S∗k and any w ∈W we have
|y − w| ≤ |y − x|+ |w − x| ≤ 3 · 2kρ.
Lemma 5.3 yields∑
|G(x, y)−G(x, z)| ≤ 3C12kρ
Hx(w).
By Lemma 5.5 we have
w∈W Hx(w) ≤ C22
−k, so the above sum is at most
3C1C2ρ. Since the union of shells S∗0 ,S
1 , . . . ,S
dlog2(r/ρ)e
covers all of Br ex-
cept for those points y within distance ρ of x, and
|y−x|≤ρ
z∼y |G(x, y)−
G(x, z)| ≤ C1ρ by Lemma 5.3, the result follows.
Proof of Theorem 1.1, Inner Estimate. Let κ and R be defined as in Lemma
5.1. Since the net number of particles to enter a site x 6= o is at most one, we
have 2d div κ(x) ≥ −1. Likewise 2d div κ(o) = n−1. Taking the divergence
in (15), we obtain
∆u(x) ≤ 1 + divR(x), x 6= o; (22)
∆u(o) = 1− n+ divR(o). (23)
Let T be the first exit time from Br, and define
f(x) = Exu(XT )− ExT + n Ex #{j < T |Xj = 0}.
Then ∆f(x) = 1 for x ∈ Br − {o} and ∆f(o) = 1 − n. Moreover f ≥ 0 on
∂Br. It follows from Lemma 2.2 with m = n that f ≥ γ − C4 on Br for a
suitable constant C4.
We have
u(x)− Exu(XT ) =
u(Xk∧T )− u(X(k+1)∧T )
Each summand on the right side is zero on the event {T ≤ k}, hence
u(Xk∧T )− u(X(k+1)∧T ) | Fk∧T
= −∆u(Xk)1{T>k}.
Taking expectations and using (22) and (23), we obtain
u(x)− Exu(XT ) ≥
1{T>k}(n1{Xk=o} − 1− divR(Xk))
= n Ex #{k < T |Xk = o} − ExT −
1{T>k}divR(Xk)
hence
u(x)− f(x) ≥ −
1{T>k} ∑
R(Xk, z)
 . (24)
Since random walk exits Br with probability at least 12d every time it reaches
a site adjacent to the boundary ∂Br, the expected time spent adjacent to
the boundary before time T is at most 2d. Since |R| ≤ 4d, the terms in (24)
with z ∈ ∂Br contribute at most 16d3 to the sum. Thus
u(x)− f(x) ≥ −
 ∑
y,z∈Br
1{T>k}∩{Xk=y}R(y, z)
− 8d2.
For y ∈ Br we have {Xk = y} ∩ {T > k} = {Xk∧T = y}, hence
u(x)− f(x) ≥ −
y,z∈Br
Px(Xk∧T = y)R(y, z)− 8d2. (25)
Write pk(y) = Px(Xk∧T = y). Note that since ∇u and κ are antisymmet-
ric, R is antisymmetric. Thus∑
y,z∈Br
pk(y)R(y, z) = −
y,z∈Br
pk(z)R(y, z)
y,z∈Br
pk(y)− pk(z)
R(y, z).
Summing over k and using the fact that |R| ≤ 4d, we conclude from (25)
u(x) ≥ f(x)−
y,z∈Br
|G(x, y)−G(x, z)| − 8d2,
where G = GBr is the Green’s function for simple random walk stopped on
exiting Br. By Lemma 5.6 we obtain
u(x) ≥ f(x)− C3(r + 1− |x|) log
r + 1− |x|
− 8d2.
Using the fact that f ≥ γ − C4, we obtain from Lemma 2.1
u(x) ≥ (r − |x|)2 − C3(r + 1− |x|) log
r + 1− |x|
The right side is positive provided r/3 ≤ |x| < r−C5 log r. For x ∈ Br/3, by
Lemma 2.3 we have u(x) > r2/4−C3r log 32 > 0, hence Br−C5 log r ⊂ An.
5.2 Outer Estimate
The following result is due to Holroyd and Propp (unpublished); we include
a proof for the sake of completeness. Notice that the bound in (26) does not
depend on the number of particles.
Proposition 5.7. Let Γ be a finite connected graph, and let Y ⊂ Z be sub-
sets of the vertex set of Γ. Let s be a nonnegative integer-valued function on
the vertices of Γ. Let Hw(s, Y ) be the expected number of particles stopping
in Y if s(x) particles start at each vertex x and perform independent sim-
ple random walks stopped on first hitting Z. Let Hr(s, Y ) be the number of
particles stopping in Y if s(x) particles start at each vertex x and perform
rotor-router walks stopped on first hitting Z. Let H(x) = Hw(1x, Y ). Then
|Hr(s, Y )−Hw(s, Y )| ≤
|H(u)−H(v)| (26)
independent of s and the initial positions of the rotors.
Proof. For each vertex u /∈ Z, arbitrarily choose a neighbor η(u). Or-
der the neighbors η(u) = v1, v2, . . . , vd of u so that the rotor at u points
to vi+1 immediately after pointing to vi (indices mod d). We assign weight
w(u, η(u)) = 0 to a rotor pointing from u to η(u), and weight w(u, vi) =
H(u) − H(vi) + w(u, vi−1) to a rotor pointing from u to vi. These assign-
ments are consistent since H is a harmonic function:
i(H(u)−H(vi)) = 0.
Figure 6: Diagram for the proof of Lemma 5.8.
We also assign weight H(u) to a particle located at u. The sum of rotor
and particle weights in any configuration is invariant under the operation of
routing a particle and rotating the corresponding rotor. Initially, the sum of
all particle weights is Hw(s, Y ). After all particles have stopped, the sum of
the particle weights is Hr(s, Y ). Their difference is thus at most the change
in rotor weights, which is bounded above by the sum in (26).
For ρ ∈ Z let
Sρ = {x ∈ Zd : ρ ≤ |x| < ρ+ 1}. (27)
Bρ = {x ∈ Zd : |x| < ρ} = S0 ∪ . . . ∪ Sρ−1.
Note that for simple random walk started in Bρ, the first exit time of Bρ
and first hitting time of Sρ coincide. Our next result is a modification of
Lemma 5(b) of [10].
Lemma 5.8. Fix ρ ≥ 1 and y ∈ Sρ. For x ∈ Bρ let H(x) = Px(XT = y),
where T is the first hitting time of Sρ. Then
H(x) ≤
|x− y|d−1
for a constant J depending only on d.
Proof. We induct on the distance |x − y|, assuming the result holds for all
x′ with |x′ − y| ≤ 1
|x− y|; the base case can be made trivial by choosing J
sufficiently large. By Lemma 5(b) of [10], we can choose J large enough so
that the result holds provided |y| − |x| ≥ 2−d−3|x − y|. Otherwise, let H1
be the outer half-space tangent to Sρ at the point of Sρ closest to x, and let
H2 be the inner half-space tangent to the ball S̃ of radius 12 |x− y| about y,
at the point of S̃ closest to x. By Lemma 5.4 applied to these half-spaces,
the probability that random walk started at x reaches S̃ before hitting Sρ
is at most 21−d. Writing T̃ for the first hitting time of S̃ ∪ Sρ, we have
H(x) ≤
x′∈eS
Px(XeT = x′)H(x′) ≤ 21−dJ ·
|x− y|
where we have used the inductive hypothesis to bound H(x′).
The lazy random walk in Zd stays in place with probability 1
, and moves
to each of the 2d neighbors with probability 1
. We will need the following
standard result, which can be derived e.g. from the estimates in [17], section
II.12; we include a proof for the sake of completeness.
Lemma 5.9. Given u ∼ v ∈ Zd, lazy random walks started at u and v can
be coupled with probability 1−C/R before either reaches distance R from u,
where C depends only on d.
Proof. Let i be the coordinate such that ui 6= vi. To define a step of the
coupling, choose one of the d coordinates uniformly at random. If the chosen
coordinate is different from i, let the two walks take the same lazy step so
that they still agree in this coordinate. If the chosen coordinate is i, let
one walk take a step while the other stays in place. With probability 1
walks will then be coupled. Otherwise, they are located at points u′, v′ with
|u′− v′| = 2. Moreover, P
|u−u′| ≥ R
for a constant C ′ depending
only on d. From now on, whenever coordinate i is chosen, let the two walks
take lazy steps in opposite directions.
∣∣∣xi = u′i + v′i2
be the hyperplane bisecting the segment [u′, v′]. Since the steps of one walk
are reflections in H1 of the steps of the other, the walks couple when they
hit H1. Let Q be the cube of side length R/
d+2 centered at u, and let H2
be a hyperplane defining one of the faces of Q. By Lemma 5.4 with h1 = 1
and h2 = R/4
d, the probability that one of the walks exits Q before the
walks couple is at most 2d · 5
1 + 1
≤ 40 d3/2
1 + 2
Lemma 5.10. With H defined as in Lemma 5.8, we have∑
|H(u)−H(v)| ≤ J ′ log ρ
for a constant J ′ depending only on d.
Proof. Given u ∈ Bρ and v ∼ u, by Lemma 5.9, lazy random walks started
at u and v can be coupled with probability 1 − 2C/|u − y| before either
reaches distance |u − y|/2 from u. If the walks reach this distance without
coupling, by Lemma 5.8 each has still has probability at most J/|u− y|d−1
of exiting Bρ at y. By the strong Markov property it follows that
|H(u)−H(v)| ≤
|u− y|d
Summing in spherical shells about y, we obtain
|H(u)−H(v)| ≤
d−1 2CJ
≤ J ′ log ρ.
We remark that Lemma 5.10 could also be inferred from Lemma 5.8
using [9, Thm. 1.7.1] in a ball of radius |u− y|/2 about u.
To prove the outer estimate of Theorem 1.1, we will make use of the
abelian property of rotor-router aggregation. Fix a finite set Γ ⊂ Zd con-
taining the origin. Starting with n particles at the origin, at each time step,
choose a site x ∈ Γ with more than one particle, rotate the rotor at x, and
move one particle from x to the neighbor the rotor points to. After a finite
number of such choices, each site in Γ will have at most one particle, and all
particles that exited Γ will be on the boundary ∂Γ. The abelian property
says that the final configuration of particles and the final configuration of
rotors do not depend on the choices. For a proof, see [4, Prop. 4.1].
In our application, we will fix ρ ≥ r and stop each particle in rotor-router
aggregation either when it reaches an unoccupied site or when it reaches the
spherical shell Sρ. Let Nρ be the number of particles that reach Sρ during
this process. Note that at some sites in Sρ, more than one particle may have
stopped. If we let each of these extra particles in turn continue performing
rotor-router walk, stopping either when it reaches an unoccupied site or
when it hits the larger shell Sρ+h, then by the abelian property, the number
of particles that reach Sρ+h will be Nρ+h. We will show that when h is order
r1−1/d, a constant fraction of the particles that reach Sρ find unoccupied sites
before reaching Sρ+h.
Proof of Theorem 1.1, Outer Estimate. Fix integers ρ ≥ r and h ≥ 1. In
the setting of Proposition 5.7, let Γ be the lattice ball Bρ+h+1, and let
Z = Sρ+h. Fix y ∈ Sρ+h and let Y = {y}. For x ∈ Sρ, let s(x) be the
number of particles stopped at x if each particle in rotor-router aggregation
is stopped either when it reaches an unoccupied site or when it reaches Sρ.
Write
H(x) = Px(XT = y)
where T is the first hitting time of Sρ+h. By Lemma 5.8 we have
Hw(s, y) =
s(x)H(x) ≤
where
is the number of particles that ever visit the shell Sρ.
By Lemma 5.10 the sum in (26) is at most J ′ log h, hence from Proposi-
ton 5.7 and (29) we have
Hr(s, y) ≤
+ J ′ log h. (30)
Let ρ(0) = r, and define ρ(i) inductively by
ρ(i+ 1) = min
ρ(i) +N2/(2d−1)
, min{ρ > ρ(i)|Nρ ≤ Nρ(i)/2}
. (31)
Fixing h < ρ(i+ 1)− ρ(i), we have
hd−1 log h ≤ N
logNρ(i) ≤ Nρ(i);
so (30) with ρ = ρ(i) simplifies to
Hr(s, y) ≤
CNρ(i)
where C = J + J ′.
Since all particles that visit Sρ(i)+h during rotor-router aggregation must
pass through Sρ(i), we have by the abelian property
Nρ(i)+h ≤
y∈Sρ(i)+h
Hr(s, y). (33)
Let Mk = #(An∩Sk). There are at most Mρ(i)+h nonzero terms in the sum
on the right side of (33), and each term is bounded above by (32), hence
Mρ(i)+h ≥ Nρ(i)+h
CNρ(i)
where the second inequality follows from Nρ(i)+h ≥ Nρ(i)/2. Summing over
h, we obtain
ρ(i+1)−1∑
ρ=ρ(i)+1
(ρ(i+ 1)− ρ(i)− 1)d. (34)
The left side is at most Nρ(i), hence
ρ(i+ 1)− ρ(i) ≤ (2dCNρ(i))
1/d ≤ N2/(2d−1)
providedNρ(i) ≥ C ′ := (2dC)2d−1. Thus the minimum in (31) is not attained
by its first argument. It follows that Nρ(i+1) ≤ Nρ(i)/2, hence Nρ(a log r) < C ′
for a sufficiently large constant a.
By the inner estimate, since the ball Br−c log r is entirely occupied, we
have ∑
Mρ ≤ ωdrd − ωd(r − c log r)d
≤ cdωdrd−1 log r.
Write xi = ρ(i+ 1)− ρ(i)− 1; by (34) we have
a log r∑
xdi ≤ cdωdr
d−1 log r,
By Jensen’s inequality, subject to this constraint,
xi is maximized when
all xi are equal, in which case xi ≤ C ′′r1−1/d and
ρ(a log r) = r +
xi ≤ r + C ′′r1−1/d log r. (35)
Since Nρ(a log r) < C
′ we have Nρ(a log r)+C′ = 0; that is, no particles reach
the shell Sρ(a log r)+C′ . Taking c′ = C ′ + C ′′, we obtain from (35)
An ⊂ Br(1+c′r−1/d log r).
Figure 7: Image of the rotor-router aggregate of one million particles under
the map z 7→ 1/z2. The colors represent the rotor directions. The white
disc in the center is the image of the complement of the occupied region.
6 Concluding Remarks
A number of intriguing questions remain unanswered. Although we have
shown that the asymptotic shape of the rotor-router model is a ball, the
near perfect circularity found in Figure 1 remains a mystery. In particular,
we do not know whether an analogue of Theorem 1.3 holds for the rotor-
router model, with constant error in the radius as the number of particles
grows.
Equally mysterious are the patterns in the rotor directions evident in
Figure 1. The rotor directions can be viewed as values of the odometer
function mod 2d, but our control of the odometer is not fine enough to
provide useful information about the patterns. If the rescaled occupied
region
π/nAn is viewed as a subset of the complex plane, it appears that
the monochromatic regions visible in Figure 1, in which all rotors point
in the same direction, occur near points of the form (1 + 2z)−1/2, where
z = a + bi is a Gaussian integer (i.e. a, b ∈ Z). We do not even have
a heuristic explanation for this phenomenon. Figure 7 shows the image
of A1,000,000 under the map z 7→ 1/z2; the monochromatic patches in the
transformed region occur at lattice points.
László Lovász (personal communication) has asked whether the occu-
pied region An is simply connected, i.e. whether its complement is con-
nected. While Theorem 1.1 shows that An cannot have any holes far from
the boundary, we cannot answer his question at present.
A final question is whether our methods could be adapted to internal
DLA to show that if n = ωdrd, then with high probability Br−c log r ⊂ In,
where In is the internal DLA cluster of n particles. The current best bound is
due to Lawler [11], who proves that with high probability Br−r1/3(log r)2 ⊂ In.
Acknowledgments
The authors thank Jim Propp for bringing the rotor-router model to our
attention and for many fruitful discussions. We also had useful discussions
with Oded Schramm, Scott Sheffield and Misha Sodin. We thank Wilfried
Huss for pointing out an error in an earlier draft. Yelena Shvets helped draw
some of the figures.
References
[1] P. Bak, C. Tang and K. Wiesenfeld, Self-organized criticality: an ex-
planation of the 1/f noise, Phys. Rev. Lett. 59, no. 4 (1987), 381–384.
[2] M. Cook, The tent metaphor, available at http://paradise.caltech.
edu/~cook/Warehouse/ForPropp/.
[3] J. N. Cooper and J. Spencer, Simulating a random walk with con-
stant error, Combin. Probab. Comput. 15 (2006) 815–822. http://www.
arxiv.org/abs/math.CO/0402323.
[4] P. Diaconis and W. Fulton, A growth model, a game, an algebra, La-
grange inversion, and characteristic classes, Rend. Sem. Mat. Univ. Pol.
Torino 49 (1991) no. 1, 95–119.
http://paradise.caltech.edu/~cook/Warehouse/ForPropp/
http://paradise.caltech.edu/~cook/Warehouse/ForPropp/
http://www.arxiv.org/abs/math.CO/0402323
http://www.arxiv.org/abs/math.CO/0402323
[5] A. Fey and F. Redig, Limiting shapes for deterministic centrally seeded
growth models, J. Stat. Phys. 130 (2008), no. 3, 579–597. http://
arxiv.org/abs/math.PR/0702450.
[6] Y. Fukai and K. Uchiyama, Potential kernel for two-dimensional ran-
dom walk. Ann. Probab. 24 (1996), no. 4, 1979–1992.
[7] M. Kleber, Goldbug variations, Math. Intelligencer 27 (2005), no. 1,
55–63.
[8] P. Koosis, La plus petite majorante surharmonique et son rapport avec
l’existence des fonctions entières de type exponentiel jouant le rôle de
multiplicateurs, Ann. Inst. Fourier (Grenoble) 33 (1983), fasc. 1, 67–
[9] G. Lawler, Intersections of Random Walks, Birkhäuser, 1996.
[10] G. Lawler, M. Bramson and D. Griffeath, Internal diffusion limited
aggregation, Ann. Probab. 20, no. 4 (1992), 2117–2140.
[11] G. Lawler, Subdiffusive fluctuations for internal diffusion limited aggre-
gation, Ann. Probab. 23 (1995) no. 1, 71–86.
[12] Y. Le Borgne and D. Rossin, On the identity of the sandpile group,
Discrete Math. 256 (2002) 775–790.
[13] L. Levine, The rotor-router model, Harvard University senior thesis
(2002), http://arxiv.org/abs/math/0409407.
[14] L. Levine and Y. Peres, Spherical asymptotics for the rotor-router
model in Zd, Indiana Univ. Math. J. 57 (2008), no. 1, 431–450.
http://arxiv.org/abs/math/0503251
[15] L. Levine and Y. Peres, The rotor-router shape is spherical, Math.
Intelligencer 27 (2005), no. 3, 9–11.
[16] L. Levine and Y. Peres, Scaling limits for internal aggregation models
with multiple sources. http://arxiv.org/abs/0712.3378.
[17] T. Lindvall, Lectures on the Coupling Method, Wiley, 1992.
[18] V. B. Priezzhev, D. Dhar, A. Dhar, and S. Krishnamurthy, Eulerian
walkers as a model of self-organised criticality, Phys. Rev. Lett. 77
(1996) 5079–82.
http://arxiv.org/abs/math.PR/0702450
http://arxiv.org/abs/math.PR/0702450
http://arxiv.org/abs/math/0409407
http://arxiv.org/abs/math/0503251
http://arxiv.org/abs/0712.3378
[19] J. Propp, Three lectures on quasirandomness, available at http://
faculty.uml.edu/jpropp/berkeley.html.
[20] F. Spitzer, Principles of Random Walk, Springer, 1976.
[21] K. Uchiyama, Green’s functions for random walks on ZN , Proc. London
Math. Soc. 77 (1998), no. 1, 215–240.
[22] J. Van den Heuvel, Algorithmic aspects of a chip-firing game, Combin.
Probab. Comput. 10, no. 6 (2001), 505–529.
http://faculty.uml.edu/jpropp/berkeley.html
http://faculty.uml.edu/jpropp/berkeley.html
	Introduction
	Basic Estimate
	Divisible Sandpile
	Classical Sandpile
	Rotor-Router Model
	Inner Estimate
	Outer Estimate
	Concluding Remarks

  The rotor-router model is a deterministic analogue of random walk. It can be
used to define a deterministic growth model analogous to internal DLA. We prove
that the asymptotic shape of this model is a Euclidean ball, in a sense which
is stronger than our earlier work. For the shape consisting of $n=\omega_d r^d$
sites, where $\omega_d$ is the volume of the unit ball in $\R^d$, we show that
the inradius of the set of occupied sites is at least $r-O(\log r)$, while the
outradius is at most $r+O(r^\alpha)$ for any $\alpha > 1-1/d$. For a related
model, the divisible sandpile, we show that the domain of occupied sites is a
Euclidean ball with error in the radius a constant independent of the total
mass. For the classical abelian sandpile model in two dimensions, with $n=\pi
r^2$ particles, we show that the inradius is at least $r/\sqrt{3}$, and the
outradius is at most $(r+o(r))/\sqrt{2}$. This improves on bounds of Le Borgne
and Rossin. Similar bounds apply in higher dimensions.

Introduction
Rotor-router walk is a deterministic analogue of random walk, first studied
by Priezzhev et al. [18] under the name “Eulerian walkers.” At each site
∗supported by an NSF Graduate Research Fellowship, and NSF grant DMS-0605166;
levine(at)math.berkeley.edu
†partially supported by NSF grant DMS-0605166; peres(at)stat.berkeley.edu
Keywords: abelian sandpile, asymptotic shape, discrete Laplacian, divisible sand-
pile, growth model, internal diffusion limited aggregation, rotor-router model
2000 Mathematics Subject Classifications: Primary 60G50; Secondary 35R35
Figure 1: Rotor-router aggregate of one million particles in Z2. Each site is
colored according to the direction of its rotor.
in the integer lattice Z2 is a rotor pointing north, south, east or west. A
particle starts at the origin; during each time step, the rotor at the particle’s
current location is rotated clockwise by 90 degrees, and the particle takes a
step in the direction of the newly rotated rotor. In rotor-router aggregation,
introduced by Jim Propp, we start with n particles at the origin; each par-
ticle in turn performs rotor-router walk until it reaches a site not occupied
by any other particles. Let An denote the resulting region of n occupied
sites. For example, if all rotors initially point north, the sequence will begin
A1 = {(0, 0)}, A2 = {(0, 0), (1, 0)}, A3 = {(0, 0), (1, 0), (0,−1)}. The region
A1,000,000 is pictured in Figure 1. In higher dimensions, the model can be
defined analogously by repeatedly cycling the rotors through an ordering of
the 2d cardinal directions in Zd.
Jim Propp observed from simulations in two dimensions that the regions
An are extraordinarily close to circular, and asked why this was so [7, 19].
Despite the impressive empirical evidence for circularity, the best result
known until now [14] says only that if An is rescaled to have unit volume,
the volume of the symmetric difference of An with a ball of unit volume tends
to zero as a power of n, as n ↑ ∞. The main outline of the argument is
summarized in [15]. Fey and Redig [5] also show that An contains a diamond.
In particular, these results do not rule out the possibility of “holes” in An
far from the boundary or of long tendrils extending far beyond the boundary
of the ball, provided the volume of these features is negligible compared to
Our main result is the following, which rules out the possibility of holes
far from the boundary or of long tendrils in the rotor-router shape. For
r ≥ 0 let
Br = {x ∈ Zd : |x| < r}.
Theorem 1.1. Let An be the region formed by rotor-router aggregation in
Zd starting from n particles at the origin and any initial rotor state. There
exist constants c, c′ depending only on d, such that
Br−c log r ⊂ An ⊂ Br(1+c′r−1/d log r)
where r = (n/ωd)1/d, and ωd is the volume of the unit ball in Rd.
We remark that the same result holds when the rotors evolve according
to stacks of bounded discrepancy; see the remark following Lemma 5.1.
Internal diffusion limited aggregation (“internal DLA”) is an analogous
growth model defined using random walks instead of rotor-router walks.
Starting with n particles at the origin, each particle in turn performs simple
random walk until it reaches an unoccupied site. Lawler, Bramson and
Griffeath [10] showed that for internal DLA in Zd, the occupied region An,
rescaled by a factor of n1/d, converges with probability one to a Euclidean
ball in Rd as n→∞. Lawler [11] estimated the rate of convergence. By way
of comparison with Theorem 1.1, if In is the internal DLA region formed
from n particles started at the origin, the best known bounds [11] are (up
to logarithmic factors)
Br−r1/3 ⊂ In ⊂ Br+r1/3
for all sufficiently large n, with probability one.
We also study another model which is slightly more difficult to define,
but much easier to analyze. In the divisible sandpile, each site x ∈ Zd
starts with a quantity of “mass” ν0(x) ∈ R≥0. A site topples by keeping up
Figure 2: Classical abelian sandpile aggregate of one million particles in Z2.
Colors represent the number of grains at each site.
to mass 1 for itself, and distributing the excess (if any) equally among its
neighbors. Thus if x has mass m > 1, then each of the 2d neighboring sites
gains mass (m−1)/2d when we topple x, and x is left with mass 1; if m ≤ 1,
then no mass moves when we topple x.
Note that individual topplings do not commute; however, the divisible
sandpile is “abelian” in the following sense.
Proposition 1.2. Let x1, x2, . . . ∈ Zd be a sequence with the property that
for any x ∈ Zd there are infinitely many terms xk = x. Let
uk(x) = total mass emitted by x after toppling x1, . . . , xk;
νk(x) = amount of mass present at x after toppling x1, . . . , xk.
Then uk ↑ u and νk → ν ≤ 1. Moreover, the limits u and ν are independent
of the sequence {xk}.
The abelian property can be generalized as follows: after performing
some topplings, we can add some additional mass and then continue top-
pling. The resulting limits u and ν will be the same as in the case when all
mass was initially present. For a further generalization, see [16].
The limiting function u in Proposition 1.2 is the odometer function for
the divisible sandpile. This function plays a central role in our analysis. The
limit ν represents the final mass distribution. Sites x ∈ Zd with ν(x) = 1
are called fully occupied. Proposition 1.2 is proved in section 3, along with
the following.
Theorem 1.3. For m ≥ 0 let Dm ⊂ Zd be the domain of fully occupied
sites for the divisible sandpile formed from a pile of mass m at the origin.
There exist constants c, c′ depending only on d, such that
Br−c ⊂ Dm ⊂ Br+c′ ,
where r = (m/ωd)1/d and ωd is the volume of the unit ball in Rd.
The divisible sandpile is similar to the “oil game” studied by Van den
Heuvel [22]. In the terminology of [5], it also corresponds to the h → −∞
limit of the classical abelian sandpile (defined below), that is, the abelian
sandpile started from the initial condition in which every site has a very
deep “hole.”
In the classical abelian sandpile model [1], each site in Zd has an integer
number of grains of sand; if a site has at least 2d grains, it topples, sending
one grain to each neighbor. If n grains of sand are started at the origin in
Zd, write Sn for the set of sites that are visited during the toppling process;
in particular, although a site may be empty in the final state, we include it
in Sn if it was occupied at any time during the evolution to the final state.
Until now the best known constraints on the shape of Sn in two dimen-
sions were due to Le Borgne and Rossin [12], who proved that
{x ∈ Z2 |x1 + x2 ≤
n/12− 1} ⊂ Sn ⊂ {x ∈ Z2 |x1, x2 ≤
n/2}.
Fey and Redig [5] proved analogous bounds in higher dimensions, and ex-
tended these bounds to arbitrary values of the height parameter h. This
parameter is discussed in section 4.
The methods used to prove the near-perfect circularity of the divisible
sandpile shape in Theorem 1.3 can be used to give constraints on the shape
of the classical abelian sandpile, improving on the bounds of [5] and [12].
Theorem 1.4. Let Sn be the set of sites that are visited by the classical
abelian sandpile model in Zd, starting from n particles at the origin. Write
n = ωdrd. Then for any � > 0 we have
Bc1r−c2 ⊂ Sn ⊂ Bc′1r+c′2
Figure 3: Known bounds on the shape of the classical abelian sandpile in
Z2. The inner diamond and outer square are due to Le Borgne and Rossin
[12]; the inner and outer circles are those in Theorem 1.4.
where
c1 = (2d− 1)−1/d, c′1 = (d− �)
−1/d.
The constant c2 depends only on d, while c′2 depends only on d and �.
Note that Theorem 1.4 does not settle the question of the asymptotic
shape of Sn, and indeed it is not clear from simulations whether the asymp-
totic shape in two dimensions is a disc or perhaps a polygon (Figure 2). To
our knowledge even the existence of an asymptotic shape is not known.
The rest of the paper is organized as follows. In section 2, we derive the
basic Green’s function estimates that are used in the proofs of Theorems 1.1,
1.3 and 1.4. In section 3 we prove Proposition 1.2 and Theorem 1.3 for the
divisible sandpile. In section 4 we adapt the methods of the previous section
to prove Theorem 1.4 for the classical abelian sandpile model. Section 5 is
devoted to the proof of Theorem 1.1.
2 Basic Estimate
Write (Xk)k≥0 for simple random walk in Zd, and for d ≥ 3 denote by
g(x) = Eo#{k|Xk = x}
the expected number of visits to x by simple random walk started at the
origin. This is the discrete harmonic Green’s function in Zd; it satisfies
∆g(x) = 0 for x 6= o, and ∆g(o) = −1, where ∆ is the discrete Laplacian
∆g(x) =
g(y)− g(x).
The sum is over the 2d lattice neighbors y of x. In dimension d = 2, simple
random walk is recurrent, so the expectation defining g(x) is infinite. Here
we define the potential kernel
g(x) = lim
gn(x)− gn(o) (1)
where
gn(x) = Eo#{k ≤ n|Xk = x}.
The limit defining g(x) in (1) is finite [9, 20], and g(x) has Laplacian ∆g(x) =
0 for x 6= o, and ∆g(o) = −1. Note that (1) is the negative of the usual
definition of the potential kernel; we have chosen this sign convention so
that g has the same Laplacian in dimension two as in higher dimensions.
Fix a real number m > 0 and consider the function on Zd
γ̃d(x) = |x|2 +mg(x). (2)
Let r be such that m = ωdrd, and let
γd(x) = γ̃d(x)− γ̃d(brce1) (3)
where e1 is the first standard basis vector in Zd. The function γd plays a
central role in our analysis. To see where it comes from, recall the divisible
sandpile odometer function of Proposition 1.2
u(x) = total mass emitted from x.
Let Dm ⊂ Zd be the domain of fully occupied sites for the divisible sandpile
formed from a pile of mass m at the origin. For x ∈ Dm, since each neighbor
y of x emits an equal amount of mass to each of its 2d neighbors, we have
∆u(x) =
u(y)− u(x)
= mass received by x−mass emitted by x
= 1−mδox.
Moreover, u = 0 on ∂Dm. By construction, the function γd obeys the same
Laplacian condition: ∆γd = 1 −mδo; and as we will see shortly, γd ≈ 0 on
∂Br. Since we expect the domain Dm to be close to the ball Br, we should
expect that u ≈ γd. In fact, we will first show that u is close to γd, and then
use this to conclude that Dm is close to Br.
We will use the following estimates for the Green’s function [6, 21]; see
also [9, Theorems 1.5.4 and 1.6.2].
g(x) =
log |x|+ κ+O(|x|−2), d = 2
ad|x|2−d +O(|x|−d), d ≥ 3.
Here ad =
(d−2)ωd
, where ωd is the volume of the unit ball in Rd, and κ is
a constant whose value we will not need to know. For x ∈ Zd we use |x| to
denote the Euclidean norm of x. Here and throughout the paper, constants
in error terms denoted O(·) depend only on d.
We will need an estimate for γd near the boundary of the ball Br. We
first consider dimension d = 2. From (4) we have
γ̃2(x) = φ(x)− κm+O(m|x|−2), (5)
where
φ(x) = |x|2 −
log |x|.
In the Taylor expansion of φ about |x| = r
φ(x) = φ(r)− φ′(r)(r − |x|) +
φ′′(t)(r − |x|)2 (6)
the linear term vanishes, leaving
γ2(x) =
(r − |x|)2 +O(m|x|−2) (7)
for some t between |x| and r.
In dimensions d ≥ 3, from (4) we have
γ̃d(x) = |x|2 + adm|x|2−d +O(m|x|−d).
Setting φ(x) = |x|2 +adm|x|2−d, the linear term in the Taylor expansion (6)
of φ about |x| = r again vanishes, yielding
γd(x) =
1 + (d− 1)(r/t)d
(r − |x|)2 +O(m|x|−d)
for t between |x| and r. Together with (7), this yields the following estimates
in all dimensions d ≥ 2.
Lemma 2.1. Let γd be given by (3). For all x ∈ Zd we have
γd(x) ≥ (r − |x|)2 +O
. (8)
Lemma 2.2. Let γd be given by (3). Then uniformly in r,
γd(x) = O(1), x ∈ Br+1 −Br−1.
The following lemma is useful for x near the origin, where the error term
in (8) blows up.
Lemma 2.3. Let γd be given by (3). Then for sufficiently large r, we have
γd(x) >
, ∀x ∈ Br/3.
Proof. Since γd(x) − |x|2 is superharmonic, it attains its minimum in Br/3
at a point z on the boundary. Thus for any x ∈ Br/3
γd(x)− |x|2 ≥ γd(z)− |z|2,
hence by Lemma 2.1
γd(x) ≥ (2r/3)2 − (r/3)2 +O(1) >
Lemmas 2.1 and 2.3 together imply the following.
Lemma 2.4. Let γd be given by (3). There is a constant a depending only
on d, such that γd ≥ −a everywhere.
3 Divisible Sandpile
Let µ0 be a nonnegative function on Zd with finite support. We start with
mass µ0(y) at each site y. The operation of toppling a vertex x yields the
mass distribution
Txµ0 = µ0 + α(x)∆δx
where α(x) = max(µ0(x)−1, 0) and ∆ is the discrete Laplacian on Zd. Thus
if µ0(x) ≤ 1 then Txµ0 = µ0 and no mass topples; if µ0(x) > 1 then the
mass in excess of 1 is distributed equally among the neighbors of x.
Let x1, x2, . . . ∈ Zd be a sequence with the property that for any x ∈ Zd
there are infinitely many terms xk = x. Let
µk(y) = Txk . . . Tx1µ0(y).
be the amount of mass present at y after toppling the sites x1, . . . , xk in
succession. The total mass emitted from y during this process is
uk(y) :=
j≤k:xj=y
µj−1(y)− µj(y) =
j≤k:xj=y
αj(y) (9)
where αj(y) = max(µj(y)− 1, 0).
Lemma 3.1. As k ↑ ∞ the functions uk and µk tend to limits uk ↑ u and
µk → µ. Moreover, these limits satisfy
µ = µ0 + ∆u ≤ 1.
Proof. Write M =
x µ0(x) for the total starting mass, and let B ⊂ Z
be a ball centered at the origin containing all points within L1-distance M
of the support of µ0. Note that if µk(x) > 0 and µ0(x) = 0, then x must
have received its mass from a neighbor, so µk(y) ≥ 1 for some y ∼ x. Since∑
x µk(x) = M , it follows that µk is supported on B. Let R be the radius
of B, and consider the quadratic weight
µk(x)|x|2 ≤MR2.
Since µk(xk)−µk−1(xk) = −αk(xk) and for y ∼ xk we have µk(y)−µk−1(y) =
αk(xk), we obtain
Qk −Qk−1 = αk(xk)
|y|2 − |xk|2
= αk(xk).
Summing over k we obtain from (9)
Qk = Q0 +
uk(x).
Fixing x, the sequence uk(x) is thus increasing and bounded above, hence
convergent.
Given neighboring vertices x ∼ y, since y emits an equal amount of mass
to each of its 2d neighbors, it emits mass uk(y)/2d to x up to time k. Thus
x receives a total mass of 1
y∼x uk(y) from its neighbors up to time k.
Comparing the amount of mass present at x before and after toppling, we
obtain
µk(x) = µ0(x) + ∆uk(x).
Since uk ↑ u we infer that µk → µ := µ0 + ∆u. Note that if xk = x, then
µk(x) ≤ 1. Since for each x ∈ Zd this holds for infinitely many values of k,
the limit satisfies µ ≤ 1.
A function s on Zd is superharmonic if ∆s ≤ 0. Given a function γ on
Zd the least superharmonic majorant of γ is the function
s(x) = inf{f(x) | f is superharmonic and f ≥ γ}.
The study of the least superharmonic majorant is a classical topic in analysis
and PDE; see, for example, [8]. Note that if f is superharmonic and f ≥ γ
f(x) ≥
f(y) ≥
s(y).
Taking the infimum on the left side we obtain that s is superharmonic.
Lemma 3.2. The limit u in Lemma 3.1 is given by u = s+ γ, where
γ(x) = |x|2 +
g(x− y)µ0(y)
and s is the least superharmonic majorant of −γ.
Proof. By Lemma 3.1 we have
∆u = µ− µ0 ≤ 1− µ0.
Since ∆γ = 1−µ0, the difference u−γ is superharmonic. As u is nonnegative,
it follows that u − γ ≥ s. For the reverse inequality, note that s + γ − u is
superharmonic on the domain D = {x | µ(x) = 1} of fully occupied sites
and is nonnegative outside D, hence nonnegative inside D as well.
As a corollary of Lemmas 3.1 and 3.2, we obtain the abelian property of
the divisible sandpile, Proposition 1.2, which was stated in the introduction.
We now turn to the case of a point source mass m started at the origin:
µ0 = mδo. More general starting distributions are treated in [16], where
we identify the scaling limit of the divisible sandpile model and show that
it coincides with that of internal DLA and of the rotor-router model. In
the case of a point source of mass m, the natural question is to identify
the shape of the resulting domain Dm of fully occupied sites, i.e. sites x for
which µ(x) = 1. According to Theorem 3.3, Dm is extremely close to a ball
of volume m; in fact, the error in the radius is a constant independent of m.
As before, for r ≥ 0 we write
Br = {x ∈ Zd : |x| < r}
for the lattice ball of radius r centered at the origin.
Theorem 3.3. For m ≥ 0 let Dm ⊂ Zd be the domain of fully occupied sites
for the divisible sandpile formed from a pile of size m at the origin. There
exist constants c, c′ depending only on d, such that
Br−c ⊂ Dm ⊂ Br+c′ ,
where r = (m/ωd)1/d and ωd is the volume of the unit ball in Rd.
The idea of the proof is to use Lemma 3.2 along with the basic estimates
on γ, Lemmas 2.1 and 2.2, to obtain estimates on the odometer function
u(x) = total mass emitted from x.
We will need the following simple observation.
Lemma 3.4. For every point x ∈ Dm − {o} there is a path x = x0 ∼ x1 ∼
. . . ∼ xk = o in Dm with u(xi+1) ≥ u(xi) + 1.
Proof. If xi ∈ Dm − {o}, let xi+1 be a neighbor of xi maximizing u(xi+1).
Then xi+1 ∈ Dm and
u(xi+1) ≥
= u(xi) + ∆u(xi)
= u(xi) + 1,
where in the last step we have used the fact that xi ∈ Dm.
Proof of Theorem 3.3. We first treat the inner estimate. Let γd be given by
(3). By Lemma 3.2 the function u−γd is superharmonic, so its minimum in
the ball Br is attained on the boundary. Since u ≥ 0, we have by Lemma 2.2
u(x)− γd(x) ≥ −C, x ∈ ∂Br
for a constant C depending only on d. Hence by Lemma 2.1,
u(x) ≥ (r − |x|)2 − C ′rd/|x|d, x ∈ Br. (10)
for a constant C ′ depending only on d. It follows that there is a constant c,
depending only on d, such that u(x) > 0 whenever r/3 ≤ |x| < r − c. Thus
Br−c−Br/3 ⊂ Dm. For x ∈ Br/3, by Lemma 2.3 we have u(x) ≥ r2/4−C >
0, hence Br/3 ⊂ Dm.
For the outer estimate, note that u−γd is harmonic onDm. By Lemma 2.4
we have γd ≥ −a everywhere, where a depends only on d. Since u vanishes
on ∂Dm it follows that u − γd ≤ a on Dm. Now for any x ∈ Dm with
r − 1 < |x| ≤ r, we have by Lemma 2.2
u(x) ≤ γd(x) + a ≤ c′
for a constant c′ depending only on d. Lemma 3.4 now implies that Dm ⊂
Br+c′+1.
4 Classical Sandpile
We consider a generalization of the classical abelian sandpile, proposed by
Fey and Redig [5]. Each site in Zd begins with a “hole” of depth H. Thus,
each site absorbs the first H grains it receives, and thereafter functions
normally, toppling once for each additional 2d grains it receives. If H is
negative, we can interpret this as saying that every site starts with h = −H
grains of sand already present. Aggregation is only well-defined in the regime
h ≤ 2d− 2, since for h = 2d− 1 the addition of a single grain already causes
every site in Zd to topple infinitely often.
Let Sn,H be the set of sites that are visited if n particles start at the
origin in Zd. Fey and Redig [5, Theorem 4.7] prove that
lim sup
# (Sn,H 4BH−1/dr) = 0,
where n = ωdrd, and 4 denotes symmetric difference. The following theo-
rem strengthens this result.
Theorem 4.1. Fix an integer H ≥ 2−2d. Let Sn = Sn,H be the set of sites
that are visited by the classical abelian sandpile model in Zd, starting from
n particles at the origin, if every lattice site begins with a hole of depth H.
Write n = ωdrd. Then
Bc1r−c2 ⊂ Sn,H
where
c1 = (2d− 1 +H)−1/d
and c2 is a constant depending only on d. Moreover if H ≥ 1− d, then for
any � > 0 we have
Sn,H ⊂ Bc′1r+c′2
where
c′1 = (d− �+H)
and c′2 is independent of n but may depend on d, H and �.
Note that the ratio c1/c′1 ↑ 1 as H ↑ ∞. Thus, the classical abelian
sandpile run from an initial state in which each lattice site starts with a
deep hole yields a shape very close to a ball. Intuitively, one can think of the
classical sandpile with deep holes as approximating the divisible sandpile,
whose limiting shape is a ball by Theorem 3.3. Following this intuition,
we can adapt the proof of Theorem 3.3 to prove Theorem 4.1; just one
additional averaging trick is needed, which we explain below.
Consider the odometer function for the abelian sandpile
u(x) = total number of grains emitted from x.
Let Tn = {x|u(x) > 0} be the set of sites which topple at least once. Then
Tn ⊂ Sn ⊂ Tn ∪ ∂Tn.
In the final state, each site which has toppled retains between 0 and 2d− 1
grains, in addition to the H that it absorbed. Hence
H ≤ ∆u(x) + nδox ≤ 2d− 1 +H, x ∈ Tn. (11)
We can improve the lower bound by averaging over a small box. For x ∈ Zd
Qk(x) = {y ∈ Zd : ||x− y||∞ ≤ k}
be the box of side length 2k + 1 centered at x, and let
u(k)(x) = (2k + 1)−d
y∈Qk(x)
u(y).
Write
T (k)n = {x |Qk(x) ⊂ Tn}.
Le Borgne and Rossin [12] observe that if T is a set of sites all of which
topple, the number of grains remaining in T is at least the number of edges
internal to T : indeed, for each internal edge, the endpoint that topples last
sends the other a grain which never moves again. Since the box Qk(x) has
2dk(2k + 1)d−1 internal edges, we have
∆u(k)(x) ≥
2k + 1
d+H −
(2k + 1)d
1Qk(o)(x), x ∈ T
n . (12)
The following lemma is analogous to Lemma 3.4.
Lemma 4.2. For every point x ∈ Tn adjacent to ∂Tn there is a path x =
x0 ∼ x1 ∼ . . . ∼ xm = o in Tn with u(xi+1) ≥ u(xi) + 1.
Proof. By (11) we have
u(y) ≥ u(xi).
Since u(xi−1) < u(xi), some term u(y) in the sum above must exceed u(xi).
Let xi+1 = y.
Proof of Theorem 4.1. Let
ξ̃d(x) = (2d− 1 +H)|x|2 + ng(x),
and let
ξd(x) = ξ̃d(x)− ξ̃d(bc1rce1).
Taking m = n/(2d− 1 +H) in Lemma 2.2, we have
u(x)− ξd(x) ≥ −ξd(x) ≥ −C(2d− 1 +H), x ∈ ∂Bc1r (13)
for a constant C depending only on d. By (11), u− ξd is superharmonic, so
u − ξd ≥ −C(2d − 1 + H) in all of Bc1r. Hence by Lemma 2.1 we have for
x ∈ Bc1r
u(x) ≥ (2d− 1 +H)
(c1r − |x|)2 − C ′(c1r)d/|x|d
, (14)
where C ′ depends only on d. It follows that u is positive on Bc1r−c2−Bc1r/3
for a suitable constant c2 depending only on d. For x ∈ Bc1r/3, by Lemma 2.3
we have u(x) > (2d− 1 +H)(c21r
2/4− C) > 0. Thus Bc1r−c2 ⊂ Tn ⊂ Sn.
For the outer estimate, let
ψ̂d(x) = (d− �+H)|x|2 + ng(x).
Choose k large enough so that 2k
d ≥ d− �, and define
ψ̃d(x) = (2k + 1)
y∈Qk(x)
ψ̂d(y).
Finally, let
ψd(x) = ψ̃d(x)− ψ̃d(bc′1rce1).
By (12), u(k) − ψd is subharmonic on T
n . Taking m = n/(d − � + H)
in Lemma 2.4, there is a constant a depending only on d, such that ψd ≥
−a(d+H) everywhere. Since u(k) ≤ (2d+H)(d+1)k on ∂T (k)n it follows that
u(k) − ψd ≤ a(d + H) + (2d + H)(d+1)k on T
n . Now for any x ∈ Sn with
c′1r − 1 < |x| ≤ c
1r we have by Lemma 2.2
u(k)(x) ≤ ψd(x) + a(d+H) + (2d+H)(d+1)k ≤ c̃2
for a constant c̃2 depending only on d, H and �. Then u(x) ≤ c′2 := (2k +
1)dc̃2. Lemma 4.2 now implies that Tn ⊂ Bc′1r+c′2 , and hence
Sn ⊂ Tn ∪ ∂Tn ⊂ Bc′1r+c′2+1.
We remark that the crude bound of (2d + H)(d+1)k used in the proof
of the outer estimate can be improved to a bound of order k2H, and the
final factor of (2k + 1)d can be replaced by a constant factor independent
of k and H, using the fact that a nonnegative function on Zd with bounded
Laplacian cannot grow faster than quadratically; see [16].
5 Rotor-Router Model
Given a function f on Zd, for a directed edge (x, y) write
∇f(x, y) = f(y)− f(x).
Given a function s on directed edges in Zd, write
div s(x) =
s(x, y).
The discrete Laplacian of f is then given by
∆f(x) = div∇f =
f(y)− f(x).
5.1 Inner Estimate
Fixing n ≥ 1, consider the odometer function for rotor-router aggregation
u(x) = total number exits from x by the first n particles.
We learned the idea of using the odometer function to study the rotor-router
shape from Matt Cook [2].
Lemma 5.1. For a directed edge (x, y) in Zd, denote by κ(x, y) the net
number of crossings from x to y performed by the first n particles in rotor-
router aggregation. Then
∇u(x, y) = −2dκ(x, y) +R(x, y) (15)
for some edge function R which satisfies
|R(x, y)| ≤ 4d− 2
for all edges (x, y).
Remark. In the more general setting of rotor stacks of bounded discrepancy,
the 4d− 2 will be replaced by a different constant here.
Proof. Writing N(x, y) for the number of particles routed from x to y, we
u(x)− 2d+ 1
≤ N(x, y) ≤
u(x) + 2d− 1
hence
|∇u(x, y) + 2dκ(x, y)| = |u(y)− u(x) + 2dN(x, y)− 2dN(y, x)|
≤ 4d− 2.
In what follows, C0, C1, . . . denote constants depending only on d.
Lemma 5.2. Let Ω ⊂ Zd − {o} with 2 ≤ #Ω <∞. Then∑
|y|1−d ≤ C0 Diam(Ω).
Proof. For each positive integer k, let
Sk = {y ∈ Zd : k ≤ |y| < k + 1}.
Then ∑
|y|1−d ≤ k1−d#Sk ≤ C ′0
for a constant C ′0 depending only on d. Since Ω can intersect at most
Diam(Ω) + 1 ≤ 2 Diam(Ω) distinct sets Sk, taking C0 = 2C ′0 the proof is
complete.
Lemma 5.3. Let G = GBr be the Green’s function for simple random walk
in Zd stopped on exiting Br. For any ρ ≥ 1 and x ∈ Br,∑
|x−y|≤ρ
|G(x, y)−G(x, z)| ≤ C1ρ. (16)
Proof. Let (Xt)t≥0 denote simple random walk in Zd, and let T be the first
exit time from Br. For fixed y, the function
A(x) = g(x− y)− Exg(XT − y) (17)
has Laplacian ∆A(x) = −δxy in Br and vanishes on ∂Br, hence A(x) =
G(x, y).
Let x, y ∈ Br and z ∼ y. From (4) we have
|g(x− y)− g(x− z)| ≤
|x− y|d−1
, y, z 6= x.
Using the triangle inequality together with (17), we obtain
|G(x, y)−G(x, z)| ≤ |g(x− y)− g(x− z)|+ Ex|g(XT − y)− g(XT − z)|
|x− y|d−1
w∈∂Br
Hx(w)
|w − y|d−1
where Hx(w) = Px(XT = w).
Write D = {y ∈ Br : |x− y| ≤}. Then∑
y 6=x
z 6=x
|G(x, y)−G(x, z)| ≤ C3ρ+ C2
w∈∂Br
Hx(w)
|w − y|1−d. (18)
Figure 4: Diagram for the Proof of Lemma 5.4.
Taking Ω = w − D in Lemma 5.2, the inner sum on the right is at most
C0Diam(D) ≤ 2C0ρ, so the right side of (18) is bounded above by C1ρ for
a suitable C1.
Finally, the terms in which y or z coincides with x make a negligible
contribution to the sum in (16), since for y ∼ x ∈ Zd
|G(x, x)−G(x, y)| ≤ |g(o)−g(x−y)|+Ex|g(XT −x)−g(XT −y)| ≤ C4.
Lemma 5.4. Let H1, H2 be linear half-spaces in Zd, not necessarily parallel
to the coordinate axes. Let Ti be the first hitting time of Hi. If x /∈ H1∪H2,
Px(T1 > T2) ≤
h1 + 1
where hi is the distance from x to Hi.
Proof. If one of H1, H2 contains the other, the result is vacuous. Otherwise,
let H̃i be the half-space shifted parallel to Hci by distance 2h2 in the direction
of x, and let T̃i be the first hitting time of Hi∪H̃i. Let (Xt)t≥0 denote simple
random walk in Zd, and write Mt for the (signed) distance from Xt to the
hyperplane defining the boundary of H1, with M0 = h1. Then Mt is a
martingale with bounded increments. Since ExT̃1 < ∞, we obtain from
optional stopping
h1 = ExMeT1 ≥ 2h2 Px (XeT1 ∈ H̃1)− Px(XeT1 ∈ H1),
hence
Px (XeT1 ∈ H̃1) ≤ h1 + 12h2 . (19)
Likewise, dM2t − t is a martingale with bounded increments, giving
ExT̃1 ≤ d ExM2eT1
≤ d(2h2 + 1)2 Px (XeT1 ∈ H̃1)
≤ d(h1 + 1)(2h2 + 1)
. (20)
Let T = min(T̃1, T̃2). Denoting by Dt the distance from Xt to the
hyperplane defining the boundary of H2, the quantity
D2t + (2h2 −Dt)
is a martingale. Writing p = Px(T = T̃2) we have
dh22 = EN0 = ENT ≥ p
(2h2)
2 + (1− p)dh22 − ExT
≥ (1 + p)dh22 − ExT
hence by (20)
h1 + 1
Finally by (19)
P(T1 > T2) ≤ p+ P(XeT1 ∈ H̃1) ≤ 52
h1 + 1
Lemma 5.5. Let x ∈ Br and let ρ = r + 1− |x|. Let
S∗k = {y ∈ Br : 2
kρ < |x− y| ≤ 2k+1ρ}. (21)
Let τk be the first hitting time of S∗k , and T the first exit time from Br. Then
Px(τk < T ) ≤ C22−k.
Figure 5: Diagram for the proof of Lemma 5.5.
Proof. Let H be the outer half-space tangent to Br at the point z ∈ ∂Br
closest to x. Let Q be the cube of side length 2kρ/
d centered at x. Then
Q is disjoint from S∗k , hence
Px(τk < T ) ≤ Px(T∂Q < T ) ≤ Px(T∂Q < TH)
where T∂Q and TH are the first hitting times of ∂Q and H. Let H1, . . . ,H2d
be the outer half-spaces defining the faces of Q, so that Q = Hc1 ∩ . . .∩H
By Lemma 5.4 we have
Px(T∂Q < TH) ≤
Px(THi < TH)
dist(x,H) + 1
dist(x,Hi)
2 dist(x,Hi)
Since dist(x,H) = |x − z| ≤ ρ and dist(x,Hi) = 2k−1ρ/
d, and ρ ≥ 1,
taking C2 = 20 d3/2(1 +
d)2 completes the proof.
Lemma 5.6. Let G = GBr be the Green’s function for random walk stopped
on exiting Br. Let x ∈ Br and let ρ = r + 1− |x|. Then∑
|G(x, y)−G(x, z)| ≤ C3ρ log
Proof. Let S∗k be given by (21), and let
W = {w ∈ ∂(S∗k ∪ ∂S
k) : |w − x| < 2
be the portion of the boundary of the enlarged spherical shell S∗k ∪∂S
k lying
closer to x. Let τW be the first hitting time of W , and T the first exit time
from Br. For w ∈W let
Hx(w) = Px(XτW∧T = w).
For any y ∈ S∗k and z ∼ y, simple random walk started at x must hit W
before hitting either y or z, hence
|G(x, y)−G(x, z)| ≤
Hx(w)|G(w, y)−G(w, z)|.
For any y ∈ S∗k and any w ∈W we have
|y − w| ≤ |y − x|+ |w − x| ≤ 3 · 2kρ.
Lemma 5.3 yields∑
|G(x, y)−G(x, z)| ≤ 3C12kρ
Hx(w).
By Lemma 5.5 we have
w∈W Hx(w) ≤ C22
−k, so the above sum is at most
3C1C2ρ. Since the union of shells S∗0 ,S
1 , . . . ,S
dlog2(r/ρ)e
covers all of Br ex-
cept for those points y within distance ρ of x, and
|y−x|≤ρ
z∼y |G(x, y)−
G(x, z)| ≤ C1ρ by Lemma 5.3, the result follows.
Proof of Theorem 1.1, Inner Estimate. Let κ and R be defined as in Lemma
5.1. Since the net number of particles to enter a site x 6= o is at most one, we
have 2d div κ(x) ≥ −1. Likewise 2d div κ(o) = n−1. Taking the divergence
in (15), we obtain
∆u(x) ≤ 1 + divR(x), x 6= o; (22)
∆u(o) = 1− n+ divR(o). (23)
Let T be the first exit time from Br, and define
f(x) = Exu(XT )− ExT + n Ex #{j < T |Xj = 0}.
Then ∆f(x) = 1 for x ∈ Br − {o} and ∆f(o) = 1 − n. Moreover f ≥ 0 on
∂Br. It follows from Lemma 2.2 with m = n that f ≥ γ − C4 on Br for a
suitable constant C4.
We have
u(x)− Exu(XT ) =
u(Xk∧T )− u(X(k+1)∧T )
Each summand on the right side is zero on the event {T ≤ k}, hence
u(Xk∧T )− u(X(k+1)∧T ) | Fk∧T
= −∆u(Xk)1{T>k}.
Taking expectations and using (22) and (23), we obtain
u(x)− Exu(XT ) ≥
1{T>k}(n1{Xk=o} − 1− divR(Xk))
= n Ex #{k < T |Xk = o} − ExT −
1{T>k}divR(Xk)
hence
u(x)− f(x) ≥ −
1{T>k} ∑
R(Xk, z)
 . (24)
Since random walk exits Br with probability at least 12d every time it reaches
a site adjacent to the boundary ∂Br, the expected time spent adjacent to
the boundary before time T is at most 2d. Since |R| ≤ 4d, the terms in (24)
with z ∈ ∂Br contribute at most 16d3 to the sum. Thus
u(x)− f(x) ≥ −
 ∑
y,z∈Br
1{T>k}∩{Xk=y}R(y, z)
− 8d2.
For y ∈ Br we have {Xk = y} ∩ {T > k} = {Xk∧T = y}, hence
u(x)− f(x) ≥ −
y,z∈Br
Px(Xk∧T = y)R(y, z)− 8d2. (25)
Write pk(y) = Px(Xk∧T = y). Note that since ∇u and κ are antisymmet-
ric, R is antisymmetric. Thus∑
y,z∈Br
pk(y)R(y, z) = −
y,z∈Br
pk(z)R(y, z)
y,z∈Br
pk(y)− pk(z)
R(y, z).
Summing over k and using the fact that |R| ≤ 4d, we conclude from (25)
u(x) ≥ f(x)−
y,z∈Br
|G(x, y)−G(x, z)| − 8d2,
where G = GBr is the Green’s function for simple random walk stopped on
exiting Br. By Lemma 5.6 we obtain
u(x) ≥ f(x)− C3(r + 1− |x|) log
r + 1− |x|
− 8d2.
Using the fact that f ≥ γ − C4, we obtain from Lemma 2.1
u(x) ≥ (r − |x|)2 − C3(r + 1− |x|) log
r + 1− |x|
The right side is positive provided r/3 ≤ |x| < r−C5 log r. For x ∈ Br/3, by
Lemma 2.3 we have u(x) > r2/4−C3r log 32 > 0, hence Br−C5 log r ⊂ An.
5.2 Outer Estimate
The following result is due to Holroyd and Propp (unpublished); we include
a proof for the sake of completeness. Notice that the bound in (26) does not
depend on the number of particles.
Proposition 5.7. Let Γ be a finite connected graph, and let Y ⊂ Z be sub-
sets of the vertex set of Γ. Let s be a nonnegative integer-valued function on
the vertices of Γ. Let Hw(s, Y ) be the expected number of particles stopping
in Y if s(x) particles start at each vertex x and perform independent sim-
ple random walks stopped on first hitting Z. Let Hr(s, Y ) be the number of
particles stopping in Y if s(x) particles start at each vertex x and perform
rotor-router walks stopped on first hitting Z. Let H(x) = Hw(1x, Y ). Then
|Hr(s, Y )−Hw(s, Y )| ≤
|H(u)−H(v)| (26)
independent of s and the initial positions of the rotors.
Proof. For each vertex u /∈ Z, arbitrarily choose a neighbor η(u). Or-
der the neighbors η(u) = v1, v2, . . . , vd of u so that the rotor at u points
to vi+1 immediately after pointing to vi (indices mod d). We assign weight
w(u, η(u)) = 0 to a rotor pointing from u to η(u), and weight w(u, vi) =
H(u) − H(vi) + w(u, vi−1) to a rotor pointing from u to vi. These assign-
ments are consistent since H is a harmonic function:
i(H(u)−H(vi)) = 0.
Figure 6: Diagram for the proof of Lemma 5.8.
We also assign weight H(u) to a particle located at u. The sum of rotor
and particle weights in any configuration is invariant under the operation of
routing a particle and rotating the corresponding rotor. Initially, the sum of
all particle weights is Hw(s, Y ). After all particles have stopped, the sum of
the particle weights is Hr(s, Y ). Their difference is thus at most the change
in rotor weights, which is bounded above by the sum in (26).
For ρ ∈ Z let
Sρ = {x ∈ Zd : ρ ≤ |x| < ρ+ 1}. (27)
Bρ = {x ∈ Zd : |x| < ρ} = S0 ∪ . . . ∪ Sρ−1.
Note that for simple random walk started in Bρ, the first exit time of Bρ
and first hitting time of Sρ coincide. Our next result is a modification of
Lemma 5(b) of [10].
Lemma 5.8. Fix ρ ≥ 1 and y ∈ Sρ. For x ∈ Bρ let H(x) = Px(XT = y),
where T is the first hitting time of Sρ. Then
H(x) ≤
|x− y|d−1
for a constant J depending only on d.
Proof. We induct on the distance |x − y|, assuming the result holds for all
x′ with |x′ − y| ≤ 1
|x− y|; the base case can be made trivial by choosing J
sufficiently large. By Lemma 5(b) of [10], we can choose J large enough so
that the result holds provided |y| − |x| ≥ 2−d−3|x − y|. Otherwise, let H1
be the outer half-space tangent to Sρ at the point of Sρ closest to x, and let
H2 be the inner half-space tangent to the ball S̃ of radius 12 |x− y| about y,
at the point of S̃ closest to x. By Lemma 5.4 applied to these half-spaces,
the probability that random walk started at x reaches S̃ before hitting Sρ
is at most 21−d. Writing T̃ for the first hitting time of S̃ ∪ Sρ, we have
H(x) ≤
x′∈eS
Px(XeT = x′)H(x′) ≤ 21−dJ ·
|x− y|
where we have used the inductive hypothesis to bound H(x′).
The lazy random walk in Zd stays in place with probability 1
, and moves
to each of the 2d neighbors with probability 1
. We will need the following
standard result, which can be derived e.g. from the estimates in [17], section
II.12; we include a proof for the sake of completeness.
Lemma 5.9. Given u ∼ v ∈ Zd, lazy random walks started at u and v can
be coupled with probability 1−C/R before either reaches distance R from u,
where C depends only on d.
Proof. Let i be the coordinate such that ui 6= vi. To define a step of the
coupling, choose one of the d coordinates uniformly at random. If the chosen
coordinate is different from i, let the two walks take the same lazy step so
that they still agree in this coordinate. If the chosen coordinate is i, let
one walk take a step while the other stays in place. With probability 1
walks will then be coupled. Otherwise, they are located at points u′, v′ with
|u′− v′| = 2. Moreover, P
|u−u′| ≥ R
for a constant C ′ depending
only on d. From now on, whenever coordinate i is chosen, let the two walks
take lazy steps in opposite directions.
∣∣∣xi = u′i + v′i2
be the hyperplane bisecting the segment [u′, v′]. Since the steps of one walk
are reflections in H1 of the steps of the other, the walks couple when they
hit H1. Let Q be the cube of side length R/
d+2 centered at u, and let H2
be a hyperplane defining one of the faces of Q. By Lemma 5.4 with h1 = 1
and h2 = R/4
d, the probability that one of the walks exits Q before the
walks couple is at most 2d · 5
1 + 1
≤ 40 d3/2
1 + 2
Lemma 5.10. With H defined as in Lemma 5.8, we have∑
|H(u)−H(v)| ≤ J ′ log ρ
for a constant J ′ depending only on d.
Proof. Given u ∈ Bρ and v ∼ u, by Lemma 5.9, lazy random walks started
at u and v can be coupled with probability 1 − 2C/|u − y| before either
reaches distance |u − y|/2 from u. If the walks reach this distance without
coupling, by Lemma 5.8 each has still has probability at most J/|u− y|d−1
of exiting Bρ at y. By the strong Markov property it follows that
|H(u)−H(v)| ≤
|u− y|d
Summing in spherical shells about y, we obtain
|H(u)−H(v)| ≤
d−1 2CJ
≤ J ′ log ρ.
We remark that Lemma 5.10 could also be inferred from Lemma 5.8
using [9, Thm. 1.7.1] in a ball of radius |u− y|/2 about u.
To prove the outer estimate of Theorem 1.1, we will make use of the
abelian property of rotor-router aggregation. Fix a finite set Γ ⊂ Zd con-
taining the origin. Starting with n particles at the origin, at each time step,
choose a site x ∈ Γ with more than one particle, rotate the rotor at x, and
move one particle from x to the neighbor the rotor points to. After a finite
number of such choices, each site in Γ will have at most one particle, and all
particles that exited Γ will be on the boundary ∂Γ. The abelian property
says that the final configuration of particles and the final configuration of
rotors do not depend on the choices. For a proof, see [4, Prop. 4.1].
In our application, we will fix ρ ≥ r and stop each particle in rotor-router
aggregation either when it reaches an unoccupied site or when it reaches the
spherical shell Sρ. Let Nρ be the number of particles that reach Sρ during
this process. Note that at some sites in Sρ, more than one particle may have
stopped. If we let each of these extra particles in turn continue performing
rotor-router walk, stopping either when it reaches an unoccupied site or
when it hits the larger shell Sρ+h, then by the abelian property, the number
of particles that reach Sρ+h will be Nρ+h. We will show that when h is order
r1−1/d, a constant fraction of the particles that reach Sρ find unoccupied sites
before reaching Sρ+h.
Proof of Theorem 1.1, Outer Estimate. Fix integers ρ ≥ r and h ≥ 1. In
the setting of Proposition 5.7, let Γ be the lattice ball Bρ+h+1, and let
Z = Sρ+h. Fix y ∈ Sρ+h and let Y = {y}. For x ∈ Sρ, let s(x) be the
number of particles stopped at x if each particle in rotor-router aggregation
is stopped either when it reaches an unoccupied site or when it reaches Sρ.
Write
H(x) = Px(XT = y)
where T is the first hitting time of Sρ+h. By Lemma 5.8 we have
Hw(s, y) =
s(x)H(x) ≤
where
is the number of particles that ever visit the shell Sρ.
By Lemma 5.10 the sum in (26) is at most J ′ log h, hence from Proposi-
ton 5.7 and (29) we have
Hr(s, y) ≤
+ J ′ log h. (30)
Let ρ(0) = r, and define ρ(i) inductively by
ρ(i+ 1) = min
ρ(i) +N2/(2d−1)
, min{ρ > ρ(i)|Nρ ≤ Nρ(i)/2}
. (31)
Fixing h < ρ(i+ 1)− ρ(i), we have
hd−1 log h ≤ N
logNρ(i) ≤ Nρ(i);
so (30) with ρ = ρ(i) simplifies to
Hr(s, y) ≤
CNρ(i)
where C = J + J ′.
Since all particles that visit Sρ(i)+h during rotor-router aggregation must
pass through Sρ(i), we have by the abelian property
Nρ(i)+h ≤
y∈Sρ(i)+h
Hr(s, y). (33)
Let Mk = #(An∩Sk). There are at most Mρ(i)+h nonzero terms in the sum
on the right side of (33), and each term is bounded above by (32), hence
Mρ(i)+h ≥ Nρ(i)+h
CNρ(i)
where the second inequality follows from Nρ(i)+h ≥ Nρ(i)/2. Summing over
h, we obtain
ρ(i+1)−1∑
ρ=ρ(i)+1
(ρ(i+ 1)− ρ(i)− 1)d. (34)
The left side is at most Nρ(i), hence
ρ(i+ 1)− ρ(i) ≤ (2dCNρ(i))
1/d ≤ N2/(2d−1)
providedNρ(i) ≥ C ′ := (2dC)2d−1. Thus the minimum in (31) is not attained
by its first argument. It follows that Nρ(i+1) ≤ Nρ(i)/2, hence Nρ(a log r) < C ′
for a sufficiently large constant a.
By the inner estimate, since the ball Br−c log r is entirely occupied, we
have ∑
Mρ ≤ ωdrd − ωd(r − c log r)d
≤ cdωdrd−1 log r.
Write xi = ρ(i+ 1)− ρ(i)− 1; by (34) we have
a log r∑
xdi ≤ cdωdr
d−1 log r,
By Jensen’s inequality, subject to this constraint,
xi is maximized when
all xi are equal, in which case xi ≤ C ′′r1−1/d and
ρ(a log r) = r +
xi ≤ r + C ′′r1−1/d log r. (35)
Since Nρ(a log r) < C
′ we have Nρ(a log r)+C′ = 0; that is, no particles reach
the shell Sρ(a log r)+C′ . Taking c′ = C ′ + C ′′, we obtain from (35)
An ⊂ Br(1+c′r−1/d log r).
Figure 7: Image of the rotor-router aggregate of one million particles under
the map z 7→ 1/z2. The colors represent the rotor directions. The white
disc in the center is the image of the complement of the occupied region.
6 Concluding Remarks
A number of intriguing questions remain unanswered. Although we have
shown that the asymptotic shape of the rotor-router model is a ball, the
near perfect circularity found in Figure 1 remains a mystery. In particular,
we do not know whether an analogue of Theorem 1.3 holds for the rotor-
router model, with constant error in the radius as the number of particles
grows.
Equally mysterious are the patterns in the rotor directions evident in
Figure 1. The rotor directions can be viewed as values of the odometer
function mod 2d, but our control of the odometer is not fine enough to
provide useful information about the patterns. If the rescaled occupied
region
π/nAn is viewed as a subset of the complex plane, it appears that
the monochromatic regions visible in Figure 1, in which all rotors point
in the same direction, occur near points of the form (1 + 2z)−1/2, where
z = a + bi is a Gaussian integer (i.e. a, b ∈ Z). We do not even have
a heuristic explanation for this phenomenon. Figure 7 shows the image
of A1,000,000 under the map z 7→ 1/z2; the monochromatic patches in the
transformed region occur at lattice points.
László Lovász (personal communication) has asked whether the occu-
pied region An is simply connected, i.e. whether its complement is con-
nected. While Theorem 1.1 shows that An cannot have any holes far from
the boundary, we cannot answer his question at present.
A final question is whether our methods could be adapted to internal
DLA to show that if n = ωdrd, then with high probability Br−c log r ⊂ In,
where In is the internal DLA cluster of n particles. The current best bound is
due to Lawler [11], who proves that with high probability Br−r1/3(log r)2 ⊂ In.
Acknowledgments
The authors thank Jim Propp for bringing the rotor-router model to our
attention and for many fruitful discussions. We also had useful discussions
with Oded Schramm, Scott Sheffield and Misha Sodin. We thank Wilfried
Huss for pointing out an error in an earlier draft. Yelena Shvets helped draw
some of the figures.
References
[1] P. Bak, C. Tang and K. Wiesenfeld, Self-organized criticality: an ex-
planation of the 1/f noise, Phys. Rev. Lett. 59, no. 4 (1987), 381–384.
[2] M. Cook, The tent metaphor, available at http://paradise.caltech.
edu/~cook/Warehouse/ForPropp/.
[3] J. N. Cooper and J. Spencer, Simulating a random walk with con-
stant error, Combin. Probab. Comput. 15 (2006) 815–822. http://www.
arxiv.org/abs/math.CO/0402323.
[4] P. Diaconis and W. Fulton, A growth model, a game, an algebra, La-
grange inversion, and characteristic classes, Rend. Sem. Mat. Univ. Pol.
Torino 49 (1991) no. 1, 95–119.
http://paradise.caltech.edu/~cook/Warehouse/ForPropp/
http://paradise.caltech.edu/~cook/Warehouse/ForPropp/
http://www.arxiv.org/abs/math.CO/0402323
http://www.arxiv.org/abs/math.CO/0402323
[5] A. Fey and F. Redig, Limiting shapes for deterministic centrally seeded
growth models, J. Stat. Phys. 130 (2008), no. 3, 579–597. http://
arxiv.org/abs/math.PR/0702450.
[6] Y. Fukai and K. Uchiyama, Potential kernel for two-dimensional ran-
dom walk. Ann. Probab. 24 (1996), no. 4, 1979–1992.
[7] M. Kleber, Goldbug variations, Math. Intelligencer 27 (2005), no. 1,
55–63.
[8] P. Koosis, La plus petite majorante surharmonique et son rapport avec
l’existence des fonctions entières de type exponentiel jouant le rôle de
multiplicateurs, Ann. Inst. Fourier (Grenoble) 33 (1983), fasc. 1, 67–
[9] G. Lawler, Intersections of Random Walks, Birkhäuser, 1996.
[10] G. Lawler, M. Bramson and D. Griffeath, Internal diffusion limited
aggregation, Ann. Probab. 20, no. 4 (1992), 2117–2140.
[11] G. Lawler, Subdiffusive fluctuations for internal diffusion limited aggre-
gation, Ann. Probab. 23 (1995) no. 1, 71–86.
[12] Y. Le Borgne and D. Rossin, On the identity of the sandpile group,
Discrete Math. 256 (2002) 775–790.
[13] L. Levine, The rotor-router model, Harvard University senior thesis
(2002), http://arxiv.org/abs/math/0409407.
[14] L. Levine and Y. Peres, Spherical asymptotics for the rotor-router
model in Zd, Indiana Univ. Math. J. 57 (2008), no. 1, 431–450.
http://arxiv.org/abs/math/0503251
[15] L. Levine and Y. Peres, The rotor-router shape is spherical, Math.
Intelligencer 27 (2005), no. 3, 9–11.
[16] L. Levine and Y. Peres, Scaling limits for internal aggregation models
with multiple sources. http://arxiv.org/abs/0712.3378.
[17] T. Lindvall, Lectures on the Coupling Method, Wiley, 1992.
[18] V. B. Priezzhev, D. Dhar, A. Dhar, and S. Krishnamurthy, Eulerian
walkers as a model of self-organised criticality, Phys. Rev. Lett. 77
(1996) 5079–82.
http://arxiv.org/abs/math.PR/0702450
http://arxiv.org/abs/math.PR/0702450
http://arxiv.org/abs/math/0409407
http://arxiv.org/abs/math/0503251
http://arxiv.org/abs/0712.3378
[19] J. Propp, Three lectures on quasirandomness, available at http://
faculty.uml.edu/jpropp/berkeley.html.
[20] F. Spitzer, Principles of Random Walk, Springer, 1976.
[21] K. Uchiyama, Green’s functions for random walks on ZN , Proc. London
Math. Soc. 77 (1998), no. 1, 215–240.
[22] J. Van den Heuvel, Algorithmic aspects of a chip-firing game, Combin.
Probab. Comput. 10, no. 6 (2001), 505–529.
http://faculty.uml.edu/jpropp/berkeley.html
http://faculty.uml.edu/jpropp/berkeley.html
	Introduction
	Basic Estimate
	Divisible Sandpile
	Classical Sandpile
	Rotor-Router Model
	Inner Estimate
	Outer Estimate
	Concluding Remarks
