\section{Introduction}

The sandpile model was
introduced by Bak, Tang, and Wiesenfeld (BTW)~\cite{BAK_1}
as a simple example of a slowly driven dissipative system
exhibiting self-organized criticality (SOC). Although today many systems
with SOC are known, it is considered as the prototype of such
models, and there is a huge literature devoted to it.
Its theoretical understanding is crucially related to
Dhar's discovery of its Abelian structure~\cite{DHAR_2} which
allows exact calculation of many of its properties~\cite{PRIEZ_1,IVASH_1}.
However, a complete analytical determination of the scaling behavior
of avalanches is still lacking.
Several approximation schemes, including
a random walk approach~\cite{OBUKHOV_1},
diffusion-like analogy~\cite{ZHANG_1}, renormalization
group~\cite{PIETRO_3,IVASH_2} and a graph
theory method~\cite{PRIEZ_2}
were proposed, but led to different results.
In addition, computer simulations -- which first had suggested
simple scaling behavior together with standard finite-size scaling (FSS) --
provide increasing evidence that the avalanche statistics is much
more complicated. While most recent authors agree upon a breakdown of FFS,
the detailed interpretation of their data is highly controversial
among different groups~\cite{DEMENECH_1,TEBALDI_1,CHESSA_2,LUEB_6}.

The standard FSS ansatz implies an asymptotic form
\begin{equation}
    P_a(a,L) \sim a^{-\tau_a}\, p(a/L^{\nu_a})      \label{eq:PaL}
\end{equation}
for the distribution of the number $a$ of toppled sites in an avalanche
(in other words, its "area"), where $L$ is the size of the lattice,
$p(z)$ is a universal function, and $\tau_a$ and $\nu_a$ are
critical exponents. This ansatz implies simple scaling of all
moments of $a$, $\langle a^n\rangle \sim L^{\sigma_n}$ with
$\sigma_n=\sigma_0+n\nu_a$. Similar ansatzes should, according to this
view, hold for the number of topplings $s$ (which differs from $a$
because sites can topple more than once in an avalanche) and for the
radius and duration of avalanches.
But recent investigations~\cite{DEMENECH_1,TEBALDI_1,CHESSA_2}
show that the two-dimensional BTW model may be characterized
by a multifractal behavior where different moments of $a$ are governed
by exponents $\sigma_n$ which are not linear in $n$ and are indeed
not related to each other for different $n$. Different reasons for this
have been proposed in~\cite{DEMENECH_1,TEBALDI_1} and in~\cite{CHESSA_2}.
Notice that multifractality of avalanches can be proven for certain
variants of the 1-$d$ sandpile model~\cite{ALI_1}.

Deviations from pure power laws had been seen
already in early simulations, but were usually interpreted as finite
size effects due to avalanches which touch the boundary of the lattice.
To illustrate that this is most likely not true, and that there is a
real problem with simple scaling, we show in Fig.~\ref{fig:ratio_Ps_Pa}
the ratio ${\cal P}_s(x,L)/{\cal P}_a(x,L)$ of the integrated
distributions ${\cal P}_s(x,L) = \int_x^\infty dx' P_s(x',L)$
for $D=2$ and for different values of $L$.
In these simulations, cylindrical boundary conditions were used (open
at $y=\pm L/2$ and periodic at $x=0, L$), and data were collected only
for avalanches starting at $y=0$. In this way we hope to have minimized
boundary effects. Also, since we do not make separate fits to $P_s(x,L)$
and $P_a(x,L)$, we have none of the uncertainties inherent in such fits.
Due to Eq.~(\ref{eq:PaL}) we would expect this ratio to scale as
$x^{\tau_a-\tau_s}$ for $x \ll \min\{L^{\nu_a}, L^{\nu_s}\} \approx
L^2$. According to analytical prediction~\cite{MAJUM_3}
and recent large scale simulations~\cite{CHESSA_2,LUEB_6,MANNA}, the difference
$\tau_a-\tau_s$ should be in the range 0.024 to 0.08.
The behavior seen in
Fig.~\ref{fig:ratio_Ps_Pa} is rather different. Although the curves
for different large $L$ perfectly superimpose in a wide range, they
are in this range not straight at all (as expected for a power law),
and their average slope in this universal range is much smaller.
Very small differences $\tau_a-\tau_s$ have been seen in several
simulations using small lattices~\cite{CHRIS_2,BIHAM}.
But it still disagrees with our data showing a lack of scaling
even for avalanches which do not reach the boundary of the lattice.
Most other variables show similar deviations from pure
power laws in $D=2$ when scrutinized closely.

In principle, one can expect that these deviations of scaling can be
explained by
assuming that the avalanche boundary advances like a pinned surface in a
random medium.
Unfortunately, this interpretation seems untenable.
As shown in~\cite{IVASH_3} (see also~\cite{GRASS_1}), avalanches
proceed in distinct {\it waves} of topplings. In each wave, any site
topples at most once. In the original version of the model waves overlap
in time, but they can be disentangled by a simple trick~\cite{IVASH_3}
so that at any time only a single wave propagates.  Therefore, if at all,
the arguments associated with pinning effects
should not apply to boundaries of avalanches but to the propagation of
wave boundaries, and they would suggest that waves show complex
behavior (notice that boundaries of successive
waves are not simply related to each other~\cite{PACZUSKI_1}).

But waves {\it do} behave simply, and do show
simple scaling behavior. This is indeed the main message of the present
paper. Our results extend analytical results derived
in~\cite{IVASH_3,DHAR_4,IVASH_4}
and large scale simulations made in~\cite{STELLA_1,PACZUSKI_1}.
While the behavior of avalanches is complex and badly understood
when typical avalanches are composed of many waves (which is the case
for $D=2$ and, to a much less degree, for $D=3$), the behavior of
single waves is simple and well understood.

\begin{figure}[t]
 \centerline{\psfig{file=fig_ratio.ps,width=8cm,angle=270}}
 \caption{Ratio ${\cal P}_s(x,L)/{\cal P}_a(x,L)$ of the integrated $s$- and
     $a$-distributions for two-dimensional sandpiles. According
     to the generally accepted FSS ansatz this should be a power law with
     exponent $\tau_a-\tau_s$ in the range $0.024$ to $0.08$ in the
     region where it is independent of $L$. The dashed line is $ x^{0.01}$.
 \label{fig:ratio_Ps_Pa}
 }
\end{figure}

In particular, we show in the present paper that FSS works for
waves of topplings. Since boundary avalanches (i.e.~avalanches which
start from an unstable boundary site) always consist
of single waves~\cite{IVASH_3}, it applies also to them.
Using the spanning
tree representation of waves, the equivalence between spanning trees
and loop-erased random walks, and rigorous estimations for
the latter, we determine the critical exponents of their
probability distributions for all dimensions $D\ge2$.
We also use the wave statistics to confirm recent
numerical~\cite{LUEB_5} and analytical~\cite{OBUKHOV_1,PRIEZ_99}
predictions for the upper critical dimension of the BTW sandpile
model. The upper and lower bounds for logarithmic
corrections to scaling for four-dimensional waves
are determined analytically and confirmed numerically.

We discuss the possibility for investigation of avalanche
distributions of the BTW model using the results obtained for
waves. One of the key characteristics in this study
is the average number $\langle n \rangle_{a}$ of waves
in avalanches of a given area $a$. Our simulations show that
in $D=2$ this number grows, most probably,
not slower than a power law with the exponent 0.17 (the first
conjecture of~\cite{MAJUM_3} was 1/6, and an analytical prediction
of~\cite{PRIEZ_2} is 1/4). This means that multiple
topplings substantially change the scaling behavior of avalanches
in comparison to waves, which was indicated in many previous studies
of the two-dimensional BTW model.

In $D=3$ the fraction of avalanches containing more than one wave
is much less than in two dimensions.
All accurate numerical estimations of the exponent~$\tau_a$
lead to $\tau_a=1.33$ which coincides
with the exact exponent $4/3$ for the wave distribution,
which could mean that the averaged number of waves in an avalanche in
$D=3$ grows not faster than logarithmically.
On the other hand, considering the numerical estimation
of this number we cannot exclude its slow
polynomial growth with the avalanche size.
Then, the scaling behavior of avalanches could be corrected for the
multiple topplings in large events, similar to the case $D=2$.

Finally, for $D\ge4$ the upper logarithmic bound for the averaged
number of waves in an avalanche~\cite{PRIEZ_99} implies that the
distributions of avalanches obey asymptotic behavior with the same
exponents as for waves.

This paper is organized as follows. In section~\ref{sec:btw_model}
we remind basic definitions of the BTW model.
Section~\ref{sec:waves_trees} is devoted to detailed explanation of
the construction of waves and their spanning tree representation.
In section~\ref{sec:green_function} we derive
the critical exponents of wave distributions.
In section~\ref{sec:lerw} we discuss analytical results for
the dynamical exponent and fractal dimension of waves.
Results of computer simulations are presented in
section~\ref{sec:numerics}.

\section{The BTW model}
\label{sec:btw_model}

We consider the $D$-dimensional BTW model on a hypercubic lattice
of linear size $L$ in which integer variables $z_i\ge 0$ represent
local energies. One perturbs the system by adding particles at
randomly chosen sites according to
\begin{equation}
    z_i \, \mapsto \, z_i+1\, .             \label{eq:perturbation}
\end{equation}
A site is called unstable if the corresponding energy $z_i$
exceeds the critical value $2D$. An unstable site relaxes, its
energy is decreased by $2D$ and the energy of the $2D$ nearest
neighbors (nn) is increased by one:
\begin{equation}
    z_{i}\;\to\;z_i\,-\,2D                  \label{eq:relaxation_1}
\end{equation}
\begin{equation}
    z_{\rm nn}\;\to\;z_{\rm nn}\;+\;1.      \label{eq:relaxation_2}
\end{equation}
In this way, the neighboring sites may be activated and
an avalanche of relaxations may proceed.
If a boundary site topples, one or more particles
leave the system. The avalanche of relaxations stops when
all sites are stable again.

One can introduce in a natural way different kinds of sub-avalanches,
e.g. clusters of sites toppled not less than a given number
of times~\cite{GRASS_1} or waves of topplings~\cite{IVASH_3}.
A relaxation event (an avalanche or sub-avalanche) is characterized
by its size $s$ (total number of topplings),
area $a$ (number of distinct toppled sites),
duration $t$ (number of parallel update steps until stable
configuration is reached),
and its radius $r$ (e.g. the radius of gyration or the
maximal distance between the origin and a toppled site).
The basic hypothesis of Bak {\it et al.}~\cite{BAK_1} claimed that
in the self-organized critical state the probability
distributions of values $s,a,t,r$ exhibit power-law behavior for
intermediate values of $x$,
\begin{equation}
    P_x(x) \sim x^{-\tau_x},                  \label{eq:prob_dist}
\end{equation}
with $x \in \{s,a,t,r\}$.

As we have seen, this hypothesis might not be true for complete
avalanches, but as we shall see it does hold for waves.
Scaling relations for the exponents $\tau_s, \tau_a, \tau_t$, and
$\tau_r$ can be obtained if one assumes that size, area,
duration and radius of ``typical" events scale as
powers of each other, for instance
\begin{equation}
    t \, \sim \, r^{\gamma_{tr}}.                \label{eq:gam_tr}
\end{equation}
Then the transformation law of probability distributions
$P_t(t) \mbox{d}t=P_r(r) \mbox{d}r$ leads to the scaling relation
\begin{equation}
    {\gamma_{tr}}\;=\;\frac{\tau_r-1}{\tau_t-1}.     \label{eq:gam_tau_tr}
\end{equation}

Again we should warn that there is a crucial assumption underlying
these relations, namely that conditional distributions $P_x(x|y)$
are narrow, and therefore Eq.~(\ref{eq:gam_tr}) holds with small
deviations for most events. It was proposed in~\cite{CHESSA_2} that
this might not be justified in $D=2$, and this is indeed the
main source of problems of this approach. Let us ignore this for the
moment and proceed nevertheless.

The scaling exponents~$\gamma_{xx'}$ are important for the
description of the extent of avalanches and their
propagation. For instance the exponent $\gamma_{sa}$ indicates
if multiple toppling events are relevant ($\gamma_{sa}>1$)
or irrelevant ($\gamma_{sa}=1$).
The exponent $\gamma_{ar}$ relating the avalanche area to
its radius $r$~equals the fractal dimension $D_f$ of the avalanche.
Finally, the exponent $\gamma_{tr}$ is usually identified
with the dynamical exponent~$z$.

If Eqs.~(\ref{eq:prob_dist}-\ref{eq:gam_tau_tr}) are applied to
waves, one has of course $\gamma_{sa}=1$, but $\gamma_{ar}$
and $\gamma_{tr}$ are non-trivial.
Our main result states that
Eqs.~(\ref{eq:prob_dist}-\ref{eq:gam_tau_tr})
do indeed apply to waves, together with
the FSS ansatz Eq.~(\ref{eq:PaL}).

\section{waves of topplings and their spanning tree representation}
\label{sec:waves_trees}

Dhar proved~\cite{DHAR_2} that all stable configurations
can be classified as either transient or recurrent. The
former can occur only during an initial transient period, but
are irrelevant for the infinite time dynamics.
He also formulated the so-called "burning algorithm" which,
on the one hand, allows one to distinguish the recurrent states
among all stable configurations, and, on the other hand, can
be used for constructing a spanning tree representation of any
given recurrent state. According to this algorithm, which also
proceeds in discrete time steps, any site $i$ is ``burnt"
at time $t$ if its energy $z_i$ is larger
than the number of its unburnt nearest neighbors at time $t-1$.
In a stable configuration, only some of boundary sites can satisfy
this condition at the first step, they can be interpreted
as origins of ``fire". Then the ``fire" propagates if new sites
become burnable at the second step.  In this way, we burn the
sites step by step, until no more sites can be burnt. If all
sites of the lattice are burnt, the initial configuration of
energies is recurrent, otherwise it it transient.

In order to obtain the spanning tree representation of recurrent
configurations~\cite{MAJUM_3},
each burnt site $i$ is connected by a bond to one of the sites
which had ``set it afire", i.e. which had caused its burning by
burning itself. If there is more than one such site, one uses an
arbitrary but definite set of rules where to place this bond.
In addition, one introduces a new site $\square$ (``sink") and
connects it to all boundary sites. On these connections, bonds
are placed to those sites which burn at $t=1$. Then
we can imagine the entire process to start by burning the
site $\square$ at time $t=0$, and generating a rooted tree
with root at $\square$. If the state is recurrent, this tree
spans the entire lattice.

Majumdar and Dhar~\cite{MAJUM_3} also noticed that
the condition for ``toppling" of a site is essentially
the same as the condition for ``burning": At each step,
the site $i$ topples if its energy $z_i$ is larger than
the number of those of its nearest neighbors which had not
toppled in the step before. Therefore,
the burning of a recurrent state is equivalent to a toppling
process initiated ``from the boundary". It implies that
if we add one particle to every boundary site (two particles
on each corner, etc.), each site will topple exactly once
during the ensuing avalanche.

The burning algorithm gives a one-to-one correspondence
between recurrent states and spanning trees.
This allows one to calculate the total number
of recurrent configurations,
the energy probabilities and the energy-energy
correlation functions~\cite{DHAR_2,PRIEZ_1,IVASH_1,MAJUM_1}.

The spanning tree representation can be constructed also
for a certain class of unstable configurations appearing
during an avalanche.
It was shown in~\cite{IVASH_3} that avalanches in the
BTW model can be decomposed into so called `waves
of topplings'.  According to this construction,
an avalanche is considered as a superposition
of successive sub-avalanches. After
perturbing the system at a given lattice site~$i$,
one allows it to relax, but prevents
the site~$i$ temporarily from toppling a second time.
After this first `wave' all sites are again stable except,
possibly, the site $i$. If $i$ is unstable, a second
wave is initiated by toppling it again. But a possible
third toppling is again delayed until this wave is finished,
and when it finally occurs, it triggers the third wave.
The procedure is repeated until the site~$i$ is stable. Note
that if the site $i$ is on the boundary, the avalanche stops
after first relaxation and consists of only one wave. More
generally, if the distance of $i$ to the boundary is $d$,
then any avalanche starting at $i$ can have at most $d+1$
waves.

To obtain the tree representation of a configuration
following a wave which had started at site $i$, we introduce
an auxiliary BTW model on a new lattice. In this lattice
we connect $i$ to the sink $\square$.
Since the site $i$ in the new model is on the boundary,
each avalanche starting at this site will consist of a single
wave. Each avalanche in the auxiliary model corresponds
indeed to a wave of some avalanche in the original model.
Applying the burning algorithm, we can construct a spanning
tree on the auxiliary lattice representing the recurrent
state of the new model. During this burning process, some
branches of the fire will be independent of site $i$, but
one branch will first pass from $\square$ to $i$ and then
propagate further. It is this latter branch which coincides
with the last wave of topplings in the original BTW model.

\begin{figure}[t]
 \epsfxsize=8.6cm
\epsffile{tree_wave.ps}
 \caption{Spanning 2-tree representation of a wave.
          Toppled sites are marked by heavy dots.
          The origin of the wave is marked by a circle.
          The dotted lines indicate the boundary of the system.
          The latter is considered as a single additional site
          $\square$, so that all non-toppled sites form a
          single connected tree.
}
 \label{fig:tree_bulk}
\end{figure}

Removing the bond between $i$ and $\square$ we obtain two
trees, one having the root at $i$ and the second at the sink
$\square$. The first tree represents the wave
and the second one corresponds to the sites not toppled in this
wave. The tree with the root $\square$ contains information
about the configuration of stable sites not affected by the wave.
We will call this union of two trees which cover the entire
lattice, a spanning two-component tree (or, simply, a spanning 2-tree).
According to this tree representation, waves with different
configurations of untoppled sites are counted as different.
An example of a spanning 2-tree is shown in Fig.~\ref{fig:tree_bulk}.

The rigorous proof of the above construction is given in~\cite{IVASH_3}.
A similar decomposition of avalanches into "inverse avalanches"
was proposed by Dhar and Manna~\cite{DHAR_4}.

An important fact concerning the wave statistics should be noted.
Since all recurrent states of the auxiliary BTW model have equal
probability of occurrence~\cite{DHAR_2}, all waves in the original
BTW model are also equally likely.

\section{Critical exponents and Green functions}
\label{sec:green_function}

Using the graph representation of waves, we can express
their probability distribution by the lattice Green function.
Consider avalanches initiated by adding a particle
at the site $i$ and spanning 2-trees with roots at $i$ and
$\square$ representing waves of these avalanches.
It was proven in~\cite{IVASH_3} that the Green function
$G_{ij}$  is related to the number
${N}_{ (\square)(ij)}$ of spanning 2-trees
where the site $j$ is in the same component as $i$:
\begin{equation}
   G_{ij}=\frac{N_{(\square)(ij)}} {N_{(\square)}},
            \label{eq:green_two_rooted}
\end{equation}
where ${N}_{ (\square)}$ is the total number of spanning trees
with the root at the site $\square$.

Consider avalanches initiated by adding a particle at the site
$i$. Since every recurrent state together with the perturbed
site $i$ completely defines the relaxation process,
the number of different possible avalanches ${N}^{(a)}_{(i)}$
started at fixed point $i$ equals to the number of recurrent
states (or the number of spanning trees ${N}_{ (\square)}$).
The number ${N}^{\scriptscriptstyle(w)}_{(ij)}$ of different waves started at
the site $i$ and covering the site $j$ corresponds to the
number ${N}_{(\square)(ij)}$ of
spanning 2-trees having sites $i,j$ on one of its components.
We can therefore rewrite Eq.~(\ref{eq:green_two_rooted}) as
\begin{equation}
    G_{ij}=\frac{N^{\scriptscriptstyle(w)}_{(ij)}}{N^{(a)}_{(i)}}.
                                  \label{eq:green_waves}
\end{equation}
Eq.~(\ref{eq:green_waves}) is another formulation of the
known result of Dhar~\cite{DHAR_2} that the expected number of
topplings at site $j$ due to adding a particle at the site $i$
is given by the Green function $G_{ij}$.

Due to uniformness of the wave statistics mentioned at
the end of section III, the probability  that a wave $W(i)$
starting at the site $i$ covers the site $j$ is equal to the
fraction of waves having this property:
\begin{equation}
     P(j \in W(i))=\frac{N^{\scriptscriptstyle(w)}_{(ij)}}{N^{\scriptscriptstyle(w)}_{(i)}},
                              \label{eq:prob_waves_01}
\end{equation}
where ${N}^{\scriptscriptstyle(w)}_{(i)}$ denotes the total number of
waves starting at $i$.

Combining Eqs.~(\ref{eq:green_waves}) and (\ref{eq:prob_waves_01}),
we write
\begin{equation}
    P(j \in W(i))=\frac{N^{\scriptscriptstyle(w)}_{(ij)}}
      {N^{(a)}_{(i)}} \frac{N^{(a)}_{(i)}}{N^{\scriptscriptstyle(w)}_{(i)}}=
      \frac{G_{ij}}{G_{ii}}=\frac{G(r)}{G(0)}.
                                        \label{eq:prob_waves_02}
\end{equation}
where we use the notation $G(r)$ for the Green function
$G_{ij}$ if the points $i$ and $j$ are separated
by the distance $r$.

On the other hand, this probability can be represented as
\begin{equation}
    P(j \in W(i)) \; = \; \int_r^\infty P^{\scriptscriptstyle(w)}(R)\;
                  \rho_R(r)\;dR\;,      \label{eq:prob_topp_02}
\end{equation}
where $P^{\scriptscriptstyle(w)}(R) $ is the
probability that the linear extent of a wave is $R$,
and $\rho_R(r)$ denotes the density of sites covered by
such a wave. The density $\rho_R(r)$ tends to 1 for large $R$
if waves are compact and isotropic, and is a function of $r$
if waves are fractal. Asymptotically for $R\gg r$, it varies as
\begin{equation}
    \rho_R(r)\approx \rho(r) \sim r^{d_f-D}.            \label{15}
\end{equation}
where $d_f$ is the fractal dimension of waves and $D$ is
the Euclidean dimension of the lattice.

Suppose that the probability distribution of the wave radius $r$
has a power law asymptotics similar to that for
avalanches [Eq.\,(\ref{eq:prob_dist})],
\begin{equation}
     P^{\scriptscriptstyle(w)}_r(r) \; \sim  \;
     r^{-\tau^{\scriptscriptstyle(w)}_r}.
                                   \label{eq:prob_wave_r}
\end{equation}
Then, the probability distribution $
P^{\scriptscriptstyle(w)}_{r}(r^\prime > r)$
scales with the exponent $\tau^{\scriptscriptstyle(w)}_r-1$.
Using Eq.\,(\ref{15})  we get the asymptotic
behavior of the probability in lhs of Eq.~(\ref{eq:prob_topp_02})
\begin{equation}
P(j \in W(i)) \; \sim \; r^{-\tau^{\scriptscriptstyle(w)}_r+1+d_f-D}.
\label{eq:green_radius}
\end{equation}

The asymptotics  of the bulk Green
function~(see for instance~\cite{JACKSON_ED}) are given by
\begin{equation}
   G(r) \; \sim \; \left \{
             \begin{array}{ll}
                        \ln{r}   & \quad {\rm for} \; D=2 \\
                                 &             \\
                        r^{2-D}  & \quad {\rm for} \; D>2 \\
             \end{array}
                   \right.
                                  \label{eq:green_asymp_bulk}
\end{equation}
which reveals that the radius exponent
$\tau^{\scriptscriptstyle(w)}_{r}$ for waves is
\begin{equation}
   \tau^{\scriptscriptstyle(w)}_r = d_f-1.
                             \label{eq:tau_radi_wave_bulk}
\end{equation}
Using Eq.\,(\ref{eq:gam_tau_tr}) we can derive
the exponents of the wave area
\begin{equation}
   \tau^{\scriptscriptstyle(w)}_a = 2 - \frac{2}{d_f}
                        \label{eq:tau_area_wave_bulk}
\end{equation}
and duration
\begin{equation}
   \tau^{\scriptscriptstyle(w)}_t = 1 + \frac{d_f-2}{z},
                        \label{eq:tau_dura_wave_bulk}
\end{equation}
respectively.

For avalanches started at a distance $b$ from the boundary,
we need the boundary Green functions which can be calculated by
the method of images:
\begin{equation}
   G({\bf r}) \; \sim \; \left \{
   \begin{array}{ll}
      \ln{|{\bf r}+{\bf b}|} -
      \ln{|{\bf r}-{\bf b}|}  & \quad {\rm for} \; D=2 \\
                              & \\
      |{\bf r}-{\bf b}|^{2-D} -
      |{\bf r}+{\bf b}|^{2-D}  & \quad {\rm for} \; D>2 \\
   \end{array}
   \right.
                     \label{eq:green_asymp_bound_01}
\end{equation}
where ${\bf b}$ is the vector perpendicular to the boundary.
On any ``equipotential" surface $G({\bf r})={\rm const}$
characterized by a length scale $\xi$ and a volume
$a\sim \xi^D$, this boundary Green function scales as
\begin{equation}
   G  \sim  b \, \xi^{1-D} \;.       \label{eq:gf_bound}
\end{equation}

If we now replace Eq.~(\ref{15}) (which is appropriate only
for isotropic cases) by its generalization $\rho \sim
\xi^{d_f-D}$, we arrive at the exponents for waves
starting near the boundary:
\begin{equation}
   \tau^{\scriptscriptstyle({\rm boundary})}_r = d_f,
               \label{eq:tau_radi_wave_bound}
\end{equation}
\begin{equation}
   \tau^{\scriptscriptstyle({\rm boundary})}_a = 2 - \frac{1}{d_f},
               \label{eq:tau_area_wave_bound}
\end{equation}
\begin{equation}
   \tau^{\scriptscriptstyle({\rm boundary})}_t = 1 + \frac{d_f-1}{z}
               \label{eq:tau_dura_wave_bound}
\end{equation}
(here and in the following we use superscript
$(w)$ for bulk {\it w}aves and ${(\rm boundary)}$ for boundary waves,
and use symbols without superscripts for avalanches).
We see that both the bulk and boundary wave exponents are
determined by the scaling exponents $d_f$ and $z$. The
dynamical exponent $z$ can be related to the fractal dimension
of the ``chemical path" on a spanning tree~\cite{MAJUM_3}
which, in turn, is equivalent to the
dimension of the loop-erased random walk (LERW)~\cite{MAJUM_2}.
As to the fractal dimension of waves $d_f$, it was proven
for all dimensions that a set of untoppled
sites, which are completely surrounded by toppled sites,
corresponds to a forbidden subconfiguration~\cite{CHRIS_2,GRASS_1}.
However, this fact does not prevent waves from being fractal.
They still could display either non-fractal or
fractal behavior depending on the dimensionality $D$.
In the next section, we will show that $d_f$
is also closely related to properties of LERW, more precisely
to the intersection probability between a LERW and a simple
random walk.

\section{Loop erased random walks, dynamical exponent and
fractal dimension of waves}
\label{sec:lerw}

In this section, we derive analytical estimates for critical
exponents of waves using their spanning tree representation and
equivalence between a chemical path on a spanning tree and
LERW.

Consider an unrestricted  random walk  on a hypercubic lattice.
The LERW introduced by Lawler~\cite{LAWLER_1} is obtained from
the simple random walk by deleting all loops along the path.
The chemical path between two sites of a tree
is the unique path along the tree edges connecting these
sites. Majumdar~\cite{MAJUM_2} has shown that the chemical path
on a spanning tree is statistically equivalent  to the LERW,
i.e., the average distance $r$ between the starting point and the
position after $l$ steps scales as $r \sim l^{\nu}$
with the same exponent $\nu$ for both of them.

In $D=2$, the exponent $\nu=4/5$ is known exactly from conformal
field theory~\cite{CONIGLIO_1}. In $D=3$, numerical estimates yield
$\nu\approx0.616$~\cite{GUTTMANN_1,BRADLEY_1}.
In $D = 4$ which is the upper critical dimension
for the LERW, $\nu =1/2$ and the simple scaling
law has logarithmic corrections~\cite{LAWLER_BOOK}.
For $D > 4$, the scaling behavior of the LERW
and chemical paths is given by the trivial value~$\nu=1/2$,
as the effects of self-intersections become
negligible above the upper critical dimension.

Returning to the BTW model, we notice that sites which topple at a
given step of wave propagation coincide with sites deleted at the
same step of the burning process, if it is started at the origin of
the wave. Since the burning process generates a tree,
there exist a unique path from the root~$i$ of the
tree (the site where the wave is initiated)
to one of the last toppled sites $i_f$  of the wave.
The number of update steps is given by the number of edges in
this path. Thus, the duration~$t$ of the
wave equals to $l$, the length of the chemical path from~$i$
to~$i_f$ on the tree and, therefore,
the dynamical exponent $\gamma_{tr}$ of waves is given by
$z=1/\nu$.

In order to find the fractal dimension of waves $d_f$, we use the
proposition proved in~\cite{PRIEZ_99}.
For this we take a site $k$ at distance $r=|k-i|$ from $i$ and
a site $j$ at distance $R=|j-i|>r$, together with some paths
$\Gamma(i,j)$ connecting $i$ with $j$ and $\Gamma(k,\square)$
connecting $k$ with the sink $\square$ (see
Fig.~\ref{fig:escape_prob}). Then the density of sites at
distance $r$ from $i$, in waves of radius $R>r$ starting at
site $i$, is given by~\cite{PRIEZ_99}:
\begin{equation}
   \rho_R(r) = \overline{P_{\rm int}(k \square|i j)},
                           \label{25}
\end{equation}
where $\overline{P_{\rm int}(k \square|ij)}$ is the probability
that $\Gamma(k,\square)$ intersects the path $\Gamma(i,j)$,
averaged over all $j$, all paths from $i$ to $j$, all $k$,
and all $\Gamma(k,\square)$.

\begin{figure}[t]
 \epsfxsize=8.6cm
 \epsffile{tree_rho_01.ps}
 \caption{A sketch of a wave initiated at the site~$i$
          and containing the site~$j$.
          The chemical path between the two sites on the
          tree representing the wave is~$\Gamma(i,j)$.
          The second path $\Gamma(k,\Box)$ corresponds to
          a random walk which starts at the site ~$k$
          and escapes the
          path~$\Gamma(i,j)$ until it is trapped at the sink.
         }
 \label{fig:escape_prob}
\end{figure}

Using the known estimations of the intersection
probabilities~\cite{LAWLER_BOOK} we can obtain from Eq.~(\ref{25})
the following upper bounds.
For $D>4$, we have
\begin{equation}
   \rho(r) = \lim_{R\to\infty} \rho_R(r) \leq C_1 r^{4-D},
                \label{26}
\end{equation}
From Eq.~(\ref{15}), we can see that the fractal dimension
of waves $d_f \leq 4$ for all $D \geq 4$.

For $D=4$, the upper bound reads
\begin{equation}
   \rho_R(r) \leq C_2 \frac{\ln (1+\alpha)}{\ln r}\;,
               \quad \alpha={R^2\over r^2}\;,       \label{27}
\end{equation}
while a lower bound was obtained in~\cite{PRIEZ_99}:
\begin{equation}
   \rho_R(r) \geq 1 - C_3 \frac{(\ln r)^{1/2}}{(\ln R)^{1/3}}.
                                              \label{28}
\end{equation}
We can see that $\rho_R(r)$ approaches 1 when
$R \rightarrow \infty,\; r$ fixed. The only fractal dimension
which is consistent with both upper and lower bounds
Eqs.~(\ref{27},\ref{28})
is 4, but there are logarithmic corrections.

For $D<4$ the lower bound Eq.~(\ref{28}) becomes stronger because of
increasing intersection probability $\overline{P_{\rm int}(k \square|i j)}$.
Thus, we conclude that the fractal dimension of waves is
\begin{equation}
   d_f \; = \; \left \{ \begin{array}{ll}
                 D  & \quad {\rm for} \; D\leq4 \\
                    &         \\
                 4  & \quad {\rm for} \; D>4 \\
               \end{array} \right.
                                          \label{32}
\end{equation}
which means that the upper critical dimension for waves is 4.

Therefore we can calculate from
Eqs.~(\ref{eq:tau_radi_wave_bulk}-\ref{eq:tau_dura_wave_bulk},
\ref{eq:tau_radi_wave_bound}-\ref{eq:tau_dura_wave_bound})
the exact values of all exponents for bulk and boundary waves
for all dimensions $D\ge 2$, with a single exception. This
exception is the exponent of the duration distribution in $D=3$,
for which we need the value of $\nu_{\scriptscriptstyle \rm LERW}(D=3)$ which is
not known exactly.

\section{Comparison with Numerical Simulations}
\label{sec:numerics}
\subsection{D=2}
\label{subsec:two_dim}

In this subsection we present the results of numerical
simulations of bulk and boundary waves in $D=2$. For these
the standard FSS works well. For any of the observables
$x=a,t,$ and $r$ it can be written as
\begin{equation}
   P_x(x,L) \; = \; L^{-\beta_x} \, g_x ( x L^{-\nu_x} ),
                                         \label{eq:fss_simple}
\end{equation}
with $\beta_x = \tau_x \nu_x$~\cite{KADANOFF_1}.

The functions $g_x(z)$ should be universal (i.e., does not
depend on the type of lattice). But for large
values of $z$ they do depend
on the type of boundary conditions (open on all four sides
or cylindrical, i.e. open in one direction and periodic in
the other) and on the aspect ratio (square or rectangle
with sides $L_1\ne L_2$). We verified that the exponents
were independent of boundary conditions and aspect ratios.
We verified also that all results were unchanged if we
threw in the sand grains with non-uniform density,
provided this density was everywhere non-zero. The latter
was very useful since it allowed us to obtain much improved
statistics from either the boundary or the central region.

\begin{figure}[t]
 \epsfysize=6.7cm
 \epsffile{a_cfss_wave_2d.eps}
 \caption{Finite-size scaling plot of the wave area distribution
          $P^{\scriptscriptstyle(w)}_a(a)$ for bulk waves in $D=2$.
          The perfect data collapse shows that $d_f=2$ and
          $\beta_a^{\scriptscriptstyle(w)}=2$, as predicted analytically.
          The dashed line represents a power law with exponent
          $\tau^{\scriptscriptstyle(w)}_a= 1$. The factor $\ln L$ comes
          from the normalization of the distributions.
          The inset verifies the scaling
          ansatz Eq.~(\protect\ref{eq:PaL}).
         }
 \label{fig:a_cfss_wave_2d}
\end{figure}

\begin{figure}[b]
 \epsfysize=6.7cm
 \epsffile{t_cfss_wave_2d.eps}
 \caption{Finite-size scaling plot of the wave duration distribution
          $P^{\scriptscriptstyle(w)}_t(t)$ for bulk waves in $D=2$.
          Here, the data collapse is achieved with
          $\beta^{\scriptscriptstyle(w)}_t=\nu^{\scriptscriptstyle(w)}_t=
          5/4$ which confirms
          again Eq.~(\protect\ref{eq:beta_estimate}).
          The dashed line represents a power law with exponent
          $\tau^{\scriptscriptstyle(w)}_t=1$. Again, the normalization factor
          $\ln L$ is needed for a good data collapse.
         }
 \label{fig:t_cfss_wave_2d}
\end{figure}

Since $d_f=2$ for $D=2$, one has $\nu_a = d_f=2$ and $\nu_r=1$,
and therefore $\nu_x=\gamma_{xr}$~\cite{LUEB_4}. The results
of the previous section, together with
$\nu_t=z=1/\nu_{\scriptscriptstyle \rm LERW}=5/4$, give
\begin{equation}
  \begin{array}{ll}
  \beta^{\scriptscriptstyle(w)}_a \; = \; 2\tau^{\scriptscriptstyle(w)}_a & = 2, \\
                   & \\
  \beta^{\scriptscriptstyle(w)}_t \; = \; z\tau^{\scriptscriptstyle(w)}_t & = 5/4, \\
  \end{array}
  \label{eq:beta_estimate}
\end{equation}
for bulk waves and
\begin{equation}
  \begin{array}{ll}
  \beta^{\scriptscriptstyle({\rm boundary})}_a \; = \; 2\tau^{\scriptscriptstyle({\rm boundary})}_a & = 3, \\
                   & \\
  \beta^{\scriptscriptstyle({\rm boundary})}_t \; = \; z\tau^{\scriptscriptstyle({\rm boundary})}_t & = 9/4, \\
  \end{array}
  \label{eq:beta_estimate1}
\end{equation}
for boundary waves.

\begin{figure}[t]
 \epsfysize=6.7cm
 \epsffile{a_cfss_wave_2d_bound.eps}
 \caption{Same as Fig.~\ref{fig:a_cfss_wave_2d}, but for boundary
          waves in $D=2$.
          Again the data collapse is obtained with
          the predicted values $\beta^{\scriptscriptstyle({\rm boundary})}_a=3$
          and d$_f=2$.
          The dashed line corresponds to
          $\tau^{\scriptscriptstyle({\rm boundary})}_a=3/2$,
          as predicted theoretically.
         }
 \label{fig:a_cfss_wave_2d_bound}
\end{figure}

\begin{figure}[b]
 \epsfysize=6.7cm
 \epsffile{t_cfss_wave_2d_bound.eps}
 \caption{Same as Fig.~\ref{fig:t_cfss_wave_2d}, but for boundary
          waves in $D=2$.
          Again the data collapse is obtained with
          the predicted values $\beta^{\scriptscriptstyle({\rm boundary})}_t=9/4$ and
          $\nu^{\scriptscriptstyle(w)}_t=5/4$.
          The dashed line corresponds to
          $\tau^{\scriptscriptstyle({\rm boundary})}_t=9/5$,
          as predicted theoretically.
         }
 \label{fig:t_cfss_wave_2d_bound}
\end{figure}

The finite-size scaling plot for the area distribution
$P^{\scriptscriptstyle(w)}_a(a)$ of bulk waves in $D=2$
is shown on Fig.~\ref{fig:a_cfss_wave_2d}.
In the inset of this Figure, as well as in the insets
of plots for other distributions of waves, we
show the collapses according to ansatz of Eq.~(\ref{eq:PaL}).
Taking $\beta^{\scriptscriptstyle(w)}_a=2$ and $d_f=2$ we see
a perfect data collapse.
The finite-size scaling ansatz of the duration
distribution is plotted in Fig.~\ref{fig:t_cfss_wave_2d}.
These data confirm that waves
are not fractal, and that their duration is as predicted by
the correspondence with spanning trees and loop-erased random
walks.

As was mentioned above, the avalanches started at the boundary
consist of a single wave, so for this type of avalanches
the distribution of avalanches coincides with that of waves.
The finite-size scaling plots for the
area and duration distributions for boundary waves
are shown in Fig.~\ref{fig:a_cfss_wave_2d_bound} and
Fig.~\ref{fig:t_cfss_wave_2d_bound}, respectively.
Again all predictions are verified, in particular we see that
$d_f=2$, i.e. also boundary waves are not fractal in $D=2$.

Finally, let us consider avalanches starting at a finite
distance from the boundary. The crossover from the boundary to
bulk behavior of the wave distribution should be described by
Eq.~(\ref{eq:green_asymp_bound_01}). More precisely, the
equipotential surfaces in $D=2$ are circles~\cite{JACKSON_ED}
with radius $\xi$
and $G \sim \ln [(b^2/\xi^2+1)^{1/2}-b/\xi]$. In the scaling
region where $\rho(\xi)=1$, we have therefore
$P^{\scriptscriptstyle(w)}_a(a|b) = -(da/d\xi)^{-1} dG/d\xi$ which gives
\begin{equation}
   P^{\scriptscriptstyle(w)}_a(a|b) \sim \; {b\over a\sqrt{a+\pi b^2}} \sim \left \{
             \begin{array}{ll}
                     b/a^{3/2}   & \; {\rm for} \; a > \pi b^2 \\
                                 &             \\
                     1/a         & \; {\rm for} \; a < \pi b^2 \\
             \end{array}
                   \right.
                                  \label{eq:Pa-boundary-2d}
\end{equation}
To check this, we simulated the BTW model with cylindrical boundary
conditions, and collected data for waves started at distance $b$
from the open boundary.
The results are plotted in Fig.~\ref{fig:cross_over_2d}.
For small $a$ we see indeed the bulk
behavior which crosses over to the boundary behavior $a^{-3/2}$
for $b^2 < a < L^2$. In the latter region we also see clearly
the linear dependence on $b$.

\begin{figure}[t]
 \epsfysize=6.6cm
 \epsffile{cross_over.eps}
 \caption{Area distributions of waves initiated
          at different distances~$b$ from the boundary
          for a system of size $256\times 1024$ with cylindrical
          boundary conditions.
          Each curve is averaged over approximately
          $10^6$ avalanches.
          The dashed lines correspond to the bulk
          ($\tau^{\scriptscriptstyle(w)}_a=1$) and boundary
          ($\tau^{\scriptscriptstyle({\rm boundary})}_a=3/2$)
          scaling behavior, respectively. The distributions
          are multiplied by $b^{-1}N(b), N(b)=\int P_a^{(w)}(a|b)da$
          in order to have the curves
          collapsing for large $a$.
          Inset: the rescaled distributions for
          $b=8,16,32,64$ demonstrate that the crossover from
          bulk to boundary behavior takes place at values of
          area of order $b^2$.
 \label{fig:cross_over_2d}}
\end{figure}

The above shows that our understanding of waves in the
two-dimensional BTW model is basically complete. In contrast, and
in spite of numerous efforts, the scaling behavior of
{\it avalanches} in the two-dimensional BTW model
is still an open problem. This is due to multiple
topplings. The average number of waves in an avalanche
scales as~\cite{DHAR_2}
\begin{equation}
   \langle n \rangle \sim \ln L \;.   \label{eq:n_w}
\end{equation}
There are also several results known
about correlations in the sizes of successive
waves~\cite{PRIEZ_2,PACZUSKI_1,KTITAREV_1,IVASH_unpub}.
Nevertheless, even most
basic questions such as the distribution of $n$ or
the dependence of $n$ on the area $a$ are not yet solved.

\begin{figure}[t]
 \epsfysize=6.1cm
 \epsffile{pn.eps}
 \caption{Probability distribution of the number of
          waves $n$ in an avalanche for $D=2$.
          The dashed line corresponds to the
          power-law behavior $P_n(n)\sim n^{-2}$ as
          predicted in~\protect\cite{MAJUM_3}.
          The data were collected from
          avalanches initiated at the center region of the
          square lattice with cylindrical boundaries.
          The inset shows the same data multiplied by $n^2$
          There, the dashed line is $\propto n^{-0.1}$.
 \label{fig:n_dist_2d}}
\end{figure}

\begin{figure}[b]
 \epsfysize=6.1cm
 \epsffile{gam_2d.eps}
 \caption{Average number of waves~$\langle n \rangle_a$
          as a function of the avalanche area~$a$ for
          $D=3$ and various system sizes $L$.
          Although the main figure looks straight at first,
          the inset displaying the rescaled average
          shows significant deviations from the
          assumed pure power-law behavior~\protect\cite{MAJUM_3}.
          The data were collected from non-dissipative
          avalanches initiated at the center region of the
          square lattice with open boundaries. Thus the curvature
          seen in the inset cannot come from avalanches reaching
          the boundary.
 \label{fig:gam_n_a_2d}}
\end{figure}

Equation (\ref{eq:n_w}) would be most easily explained if
$P_n(n)$ were simply $\sim 1/n^2$. Present data seem to agree
with this for the largest lattices (Fig.~\ref{fig:n_dist_2d}),
but actually the data are better
fited with a power $1/n^{2.1}$ than with $1/n^2$ (see
inset).
Similar results are obtained for
$\langle n\rangle_a$, the average number of waves in avalanches
with fixed $a$ (Fig.~\ref{fig:gam_n_a_2d}).
Although they seem to scale like a power of $a$,
as assumed in~\cite{MAJUM_3}, a closer study shows significant
deviations which seem hard to explain as finite size effects.

There are several recent papers~\cite{DEMENECH_1,DROSSEL_1}
which try to explain these problems by unexpected features
of avalanches which reach the boundary. But data such as
those shown in Figs.~\ref{fig:ratio_Ps_Pa},\ref{fig:gam_n_a_2d} indicate that
there are already problems with avalanches which do not
reach the boundary.

\subsection{D=3}
\label{subsec:three_dim}

Simulations of waves of the BTW model in $D=3$
also give good agreement with our analytical results
of Sections \ref{sec:green_function} and \ref{sec:lerw}.
For bulk waves we now have $\beta^{\scriptscriptstyle(w)}_a=4,\;
\tau^{\scriptscriptstyle(w)}_a = 4/3,\; \beta^{\scriptscriptstyle(w)}_t = 1/\nu_{\scriptscriptstyle \rm LERW}+1 \approx
2.623,$ and $\tau^{\scriptscriptstyle(w)}_t = 1+\nu_{\scriptscriptstyle \rm LERW} \approx 1.616$.
For boundary waves, the corresponding numbers are
$\beta^{\scriptscriptstyle({\rm boundary})}_a=5,\; \tau^{\scriptscriptstyle({\rm boundary})}_a = 5/3,\;
\beta^{\scriptscriptstyle({\rm boundary})}_t = 1/\nu_{\scriptscriptstyle \rm LERW}+2\approx 3.623$, and
$\tau^{\scriptscriptstyle({\rm boundary})}_t = 1+2\nu_{\scriptscriptstyle \rm LERW} \approx 2.232$.
For these $\beta$ values, the data collapses of bulk
(Fig.~\ref{fig:a_cfss_wave_3d},\ref{fig:t_cfss_wave_3d}) and boundary
(Fig.~\ref{fig:a_cfss_wave_3d_bound},\ref{fig:t_cfss_wave_3d_bound})
waves are perfect.
They confirm also the analytical predictions for the $\tau$ exponents,
verifying in particular that the waves have fractal dimension
$d_f=3$.

\begin{figure}[b]
 \epsfysize=6.8cm
 \epsffile{a_cfss_wave_3d.eps}
 \caption{Finite-size scaling plots of the area
          distribution $P^{\scriptscriptstyle(w)}_a(a)$ for bulk
          avalanches in $D=3$.
          Assuming compact avalanches ($d_f=3$) we get
          good data collapses and the resulting $\tau^{(\scriptscriptstyle w)}_a$
          exponent agrees with the
          theoretical prediction (dashed lines).
         }
 \label{fig:a_cfss_wave_3d}
\end{figure}

\begin{figure}[t]
 \epsfysize=6.7cm
 \epsffile{t_cfss_wave_3d.eps}
 \caption{Finite-size scaling plots of the duration
          distribution $P^{\scriptscriptstyle (w)}_t(t)$ for bulk
          avalanches in $D=3$.
          Using compact avalanches ($z=1/\nu_{\scriptscriptstyle \rm LERW}$) we get
          good data collapses and the resulting $\tau^{(\scriptscriptstyle w)}_t$
          exponent agrees with the
          theoretical prediction (dashed lines).
         }
 \label{fig:t_cfss_wave_3d}
\end{figure}

\begin{figure}[b]
 \epsfysize=6.8cm
 \epsffile{a_cfss_wave_3d_bound.eps}
\caption{Finite-size scaling plots of the area
          distribution for boundary waves in $D=3$.
          The data collapses confirm again $d_f=3$,
          and the dashed line demonstrates the agreement with
          the predicted value for $\tau_a^{({\rm boundary})}$.
         }
 \label{fig:a_cfss_wave_3d_bound}
\end{figure}

\begin{figure}[t]
 \epsfysize=6.8cm
 \epsffile{t_cfss_wave_3d_bound.eps}
\caption{Finite-size scaling plot of the duration
          distributions for boundary waves in $D=3$.
          The data collapses confirm again $z=1/\nu_{\scriptscriptstyle \rm LERW}$,
          and the dashed line demonstrates the agreement with
          the predicted value for $\tau_t^{({\rm boundary})}$.
         }
 \label{fig:t_cfss_wave_3d_bound}
\end{figure}

Due to the rarity of multiple topplings, avalanche distributions
coincide within the displayed accuracy with wave distributions,
when plotted as in Fig.~\ref{fig:a_cfss_wave_3d} and
Fig.~\ref{fig:t_cfss_wave_3d}.
In order to show significant results for multiple topplings, we have
to present the data differently. In Fig.~\ref{fig:gam_n_s_3d}
we plotted the average number of waves at fixed $a$,
$\langle n\rangle_a$, against $\log a$. Neither using a logarithmic
scale for $\langle n\rangle_a$ (main figure) nor a linear scale
(inset) give perfectly straight lines. Thus the data can
be interpreted either as a power law with a very small
exponent,
\begin{equation}
   \langle n \rangle_a \sim a^\alpha \;\;,\quad \alpha \approx 0.06
                                \label{eq:n_s_D_3}
\end{equation}
or as a logarithmic growth.

In the latter case, we would of course have
$\tau^{\scriptscriptstyle(w)}_a=\tau_a=\tau_s$. In the opposite case of
a power law with exponent $\alpha \approx 0.06$ we can give
crude estimates for the differences between these $\tau$
exponents, using some heuristic assumptions.

The first assumption is that different waves
in the same avalanche involve essentially the same sites.
If this is true, we should have
$P^{\scriptscriptstyle(w)}_a(a)da \approx \langle n \rangle_a P_a(a)da $.
Using this together with Eq.~(\ref{eq:n_s_D_3}), we
find $\tau_a = 4/3 + \alpha \approx 1.39$. Since the basic
assumption here is most likely not justified, this is
only a very crude (and most likely too large, in particular since
the growth of $\langle n \rangle_a P_a(a)da $ could be logarithmic)
estimate for
the difference between $\tau_a$ and $\tau^{\scriptscriptstyle(w)}_a$.

An estimate for the difference between $\tau_a$ and $\tau_s$
is obtained as follows. An upper bound for the size
of an avalanche of area $a$ containing $n$ waves is $s < na$.
Therefore, the assumption $s \sim \langle n \rangle_a a$ leads to the
maximal difference between the area and size exponents. Using
Eq.~(\ref{eq:n_s_D_3}) we get $a \gtrsim s^{1/(1+\alpha)}$. Then,
Eq.~(\ref{eq:gam_tau_tr}) gives $(\tau_s-1)\geq (\tau_a-1)/(1+\alpha)$
and we can conclude that
\begin{equation}
\tau_a-\tau_s \leq \frac{\alpha}{1+\alpha} (\tau_a-1) \leq 0.02.
\label{eq:estimate_tau_s}
\end{equation}

A direct verification of such slight differences between the
$\tau$ exponents could be tried by plotting ratios of the
distributions, as was done in Fig.~\ref{fig:ratio_Ps_Pa} for $D=2$.
We do not show any such ratio here, since they are all
very close to 1 in the region where the distributions
should follow power laws, and the deviations from 1 seem
not to be simple powers.
\vfill
\begin{figure}[h]
 \epsfxsize=7.2cm
 \epsffile{gam_n_a_3d.eps}
 \caption{Average number of waves~$\langle n \rangle_a$
          as a function of the avalanche area~$a$ for
          $D=3$ and various values of the system size $L$. In order to
          minimize finite size effects, cylindrical
          boundary conditions with one open and two periodic
          directions were used and one half of all sand grains
          were thrown onto the central planes $y=L/2$ and
          $y=L/2+1$. The y-axis
          is logarithmic in the main plot, and linear in the
          inset. Neither way of plotting gives perfectly
          straight lines in the region where that data for
          different $L$ collapse. Although the main figure
          looks more straight at first sight, a more careful
          inspection shows a slight downward curvature.
 \label{fig:gam_n_s_3d}}
\end{figure}

\subsection{D=4}

As the dimension of the BTW model increases, multiple
toppling events in avalanches become more and more rare.
For $D\ge 4$ it was shown in~\cite{PRIEZ_99} that
$\langle n \rangle_a$ grows not faster than logarithmically,
i.e. $\alpha=0$. As we already mentioned in the previous
subsection, this means $\tau_a^{\scriptscriptstyle(w)}=\tau_a=\tau_s$,
i.e., the scaling behavior of waves and avalanches in $D=4$ is
characterized by the same exponents and scaling functions.

At the upper critical dimension $D_{\rm u}=4$ logarithmic
corrections to the simple scaling are essential.
The probability distributions of the radius, duration, area
and the scaling relations between them are expected to have
the form (cf.~\cite{LUEB_5})
\begin{equation}
   P_r(r) \sim \frac{(\ln r)^{x_r}}{r^3}, \;
   P_t(t) \sim \frac{(\ln t)^{x_t}}{t^2}, \;
   P_a(a) \sim \frac{(\ln a)^{x_a}}{a^{3/2}},
                        \label{eq:prob_dist_4d}
\end{equation}
and
\begin{equation}
   a \; \sim \; \frac{r^4}{(\ln r)^{N_a}}, \quad
   t \; \sim \; \frac{r^2}{(\ln r)^{N_t}},
                        \label{t}
\end{equation}
respectively.

The  exponents of logarithmic corrections
$x_r$, $x_a$, $x_t$, $N_a$, and $N_t$ obey
the scaling relations~\cite{LUEB_5}
\begin{equation}
   x_r=x_a + N_a/2, \quad x_r=x_t + N_t
\label{eq:scal_rel_xr_xa_xt}
\end{equation}
which follow straightforwardly from Eqs.~(\ref{eq:prob_dist_4d})
and (\ref{t}).

Using arguments similar to those at the end of the previous subsection,
we obtain an inequality for the exponents of logarithmic corrections
for waves and avalanches
\begin{equation}
   x_a^{\scriptscriptstyle(w)} = x_a \le x_s.
\label{x}
\end{equation}

This allows us to compare below analytical estimations for waves
with numerical results for avalanches.

It follows from Lawler's results~\cite{LAWLER_BOOK} discussed in
Sect.~\ref{sec:lerw} that
$N_t=1/3$ exactly. But as with all logarithmic corrections,
a numerical verification is not easy. The main reason is that
the logarithms are never very much larger than 1, even for the
largest simulations. Therefore the next-to-leading terms
(which are typically suppressed by powers of the same
logarithms) are in general not negligible. In view of this,
the disagreement with recent simulations~\cite{LUEB_5} which
had suggested $N_t \approx 1/2$ should not be taken serious.
Data for the mean squared radius of avalanches with fixed duration $t$
are shown in Fig.~\ref{fig:gam_r_t_4d}. More precisely, since we
expect
\begin{equation}
   {\langle r\rangle_t^2 \over t} \sim (\ln t)^{1/3},
          \label{eq:lerw_scal_log}
\end{equation}
we plotted $[\langle r\rangle_t^2 /t]^3$ against $\ln t$. Apart from very
large $t$ when the finiteness of the lattice makes itself seen,
we observe essentially a straight line (which is a bit fortuitious
since there are also $1/t$ corrections which are
important for small $t$). At the same time, a power
law dependence $\langle r\rangle_t^2 \sim t^{2\nu}$ with $\nu>1/2$, as
would be expected if $D_c>4$, seems ruled out.

\begin{figure}[t]
 \epsfysize=6.8cm
  \epsffile{gam_r_t_4d.eps}
 \caption{The gyration radius of avalanches as a function of
          their duration.
          We plot the sixth power of the average rescaled gyration radius,
          $[r_g^2 t^{-1}]^3$, in a logarithmic diagram, since this should result
          in a straight line according to Eq.~(\protect\ref{eq:lerw_scal_log}).
          Such a linear regime is indeed observed (dashed line; its
          slope and intercept are not predicted by theory),
          and it increases with the system size~$L$.
 \label{fig:gam_r_t_4d}}
\end{figure}

\begin{figure}[t]
 \epsfysize=6.8cm
\epsffile{size_prediction_02.eps}
 \caption{Area distributions for avalanches in $D=4$.
          In order to render the plots more significant, first of all
          the dominant $s$-dependence was removed by multiplying
          with $s^{3/2}$. Then, in order to check whether the remaining
          $s$-dependence in the scaling region is compatible
          with logarithmic corrections as proposed in
          Eq.~(\protect\ref{eq:prob_dist_4d}), we also
          divided by powers $(\ln s)^x$ and shifted the resulting curves
          horizontally and vertically in order to avoid overlaps.
          The best agreement is found with $x\approx 0.25$.
         }
 \label{fig:a_fss_4d_log}
\end{figure}

For the other exponents $x_r, x_a, x_t, N_a$ we can only give
inequalities from the analytical results of section~\ref{sec:lerw}.
The upper bound Eq.~(\ref{27}) for
the density of waves leads to the relation
\begin{equation}
   \rho(r) \; \sim \; (\ln r)^{-\delta}, \quad \delta \ge 1.
\label{rho}
\end{equation}
Using this asymptotics and Eq.~(\ref{eq:prob_waves_02}) we get
\begin{equation}
   P_r^{\scriptscriptstyle(w)}(r) \; \sim \; G(r)/\rho(r) \;
                   \sim \; \frac{(\ln r)^\delta}{r^3},
\label{P(r)_alpha}
\end{equation}
which gives $x_r=\delta$.
The area~$a$ of a wave scales in leading order as
\begin{equation}
  a \; \sim \; \int_1^r \rho(r')\, r'^3 \, dr' \sim \frac{r^4}{(\ln r)^\delta},
\label{a_alpha}
\end{equation}
giving $N_a=\delta$.
Hence, from Eq.~(\ref{eq:scal_rel_xr_xa_xt}) we have
\begin{equation}
    x_a \; = \; \delta /2 \; \ge \; 1/2.             \label{x_a}
\end{equation}

In order to verify these predictions -- and to verify, in the
first place, that deviations from power laws with the mean
field exponents $\tau_r=3,\tau_t=2, \tau_a=3/2$ cannot be
eliminated by changing these exponents -- we performed extensive
simulations.
Numerical data of the size distribution $P_s(s)$ are shown
in Fig.~\ref{fig:a_fss_4d_log}.
Lattice sizes range from $L=32$ to $L=144$.
After multiplying with $s^{3/2}$, we see
indeed no indication that the remaining $s$-dependence follows
a clean power law. In order to find the expected logarithmic
corrections, we multiplied these data by $(\ln s)^{x_a}$, with
several trial values for the exponent $x_a$. Taken at face
value, this would give best fits with $x_a \approx 0.25$. In
view of inequality~(\ref{eq:scal_rel_xr_xa_xt}) and of the
difficulties in obtaining correct logarithmic corrections
mentioned above, we propose that indeed $x_a=1/2$. From
the relations Eq.~(\ref{eq:scal_rel_xr_xa_xt})
we get then $x_r=N_a=1$ and $x_t=2/3$.

\section{Conclusions}

We have studied probability distributions of waves of topplings
in the BTW model on $D$-dimensional rectangular lattices for
$D\ge2$. We have proved analytically that waves as well as
boundary avalanches do exhibit critical
behavior and that their probability distributions display power
law asymptotics. We have derived exact values of critical
exponents of these distributions. We have proven analytically
that the upper critical dimension of the BTW model is $D_{\rm u}=4$,
showing that previously observed deviations from mean field
behavior at $D=4$~\cite{CHRIS_2,CHESSA_1}
are due to logarithmic corrections.
All these results have been confirmed by extensive
numerical simulations. During these simulations we have also
verified that wave distributions follow the standard finite size
scaling ansatz.
The exponent of the leading logarithmic correction to the
distribution of avalanche life times (or, more precisely, lifetimes
of waves) has been derived exactly from known asymptotics of
loop-erased random walks. Estimations are given
for the exponents of the logarithmic corrections to the other
distributions.

We therefore have now a rather complete picture of the dynamics
of single waves in the BTW model for all dimensions. For $D\ge 4$
this means that we also understand avalanche dynamics, since
multiple topplings are so rare that they can be neglected. For
$D=2$ the latter is certainly not true, and our understanding
of avalanche dynamics is still incomplete.
For $D=3$, finally, multiple topplings represent a small
but not negligible effect, and we have hope that the
problem will be solved soon.

\acknowledgements
The authors thank E.\,V.~Ivashkevich for very helpful remarks.
DVK acknowledges the financial support from the Alexander von
Humboldt Foundation and the kind hospitality of the Computational
Physics group in Duisburg University,
where part of this work was done.
VBP was supported by RFBR through Grant No 99-01-00882
and by the Heisenberg-Landau program.
He also thanks NIC, Forschungszentrum J\"ulich for
the support and hospitality.