# Flowing sand - a physical realization of Directed Percolation

## 1 The Douady-Daerr experiment.

Glass beads (“sand”) of diameter $`250`$-$`425\mu `$ are poured uniformly at the top of an inclined plane (size $``$ $`1m`$), covered by a rough velvet cloth; the angle of inclination $`\phi _0`$ can be varied. As the beads flow down, a thin layer of thickness $`h=h_d(\phi _0)`$, consisting of several monolayers, settles and remains immobile. At this thickness the sand is dynamically stable; the value of $`h_d`$ decreases with increasing angle of inclination.

For each $`\phi _0`$ there exists another thickness $`h_s`$ with $`h_s(\phi _0)>h_d(\phi _0)`$, beyond which a static layer becomes unstable. Hence there exists a region in the $`(\phi ,h)`$ plane, in which a static layer is stable but a flowing one is unstable. We can now take the system, that settled at $`h_d(\phi _0)`$, and increase its angle of inclination to $`\phi =\phi _0+\mathrm{\Delta }\phi `$, staying within this region of mixed stability. The layer will not flow spontaneously, but if we disturb it, generating a flow at the top, an avalanche will propagate, leaving behind a layer of thickness $`h_d(\phi )`$. These avalanches had the shape of a fairly regular triangle with opening angle $`\theta `$. As the increment $`\mathrm{\Delta }\phi `$ decreases, the value of $`\theta `$ decreases as well, vanishing as $`\mathrm{\Delta }\phi 0`$. This calls for testing a power law behavior of the form

$$\theta (\mathrm{\Delta }\phi )^x.$$
(1)   
If instead of increasing $`\phi `$ we lower the plane, i.e., go to $`\mathrm{\Delta }\phi <0`$, the thickness of our system, $`h_d(\phi _0)`$, is less than the present thickness of dynamic stability, $`h_d(\phi )`$. In this case an initial perturbation should not propagate, it will rather die out after a certain time (or beyond a certain size $`\xi _{}`$ of the transient avalanche). As $`|\mathrm{\Delta }\phi |0`$, we expect this decay length to grow with a power law:

$$\xi _{}(\mathrm{\Delta }\phi )^\nu _{}.$$
(2)   
Hence by pouring sand at inclination $`\phi _0`$, DD produced a critical system, precisely at the borderline (with respect to changing the angle) between a stable regime $`\phi <\phi _0`$, in which perturbations die out, and an unstable one $`\phi >\phi _0`$, where perturbations persist and spread. The preparation procedure can be considered as a special kind of self-organized criticality (SOC) which differs from standard SOC models  in which a slow driving force (acting on a time scale much smaller than that of the system’s dynamic response) causes evolution to a critical state. Here avalanches are started by hand one by one.

To associate this threshold phenomenon with DP, denote by $`p`$ the percolation probability and by $`p_c`$ its critical value. We associate the change in tilt with $`pp_c`$, i.e., assume that near the angle of preparation $`\phi _0`$ the behavior of the sand system is related to a DP problem with $`\mathrm{\Delta }\phi pp_c`$. The exponent $`\nu _{}`$ should be compared with the known values for DP and CDP. The exponent $`x`$ in Eq. (1) can also be measured and compared with

$$\mathrm{tan}\theta \xi _{}/\xi _{}(\mathrm{\Delta }\phi )^{\nu _{}\nu _{}}.$$
(3)   
## 2 Definition of the model.

To write down a simple model based on the physics of the flowing sand, we adopt an observation made by DD, that in the regime of interest ($`\phi \phi _0`$) grains of the top layer rest on grains of the layers below (rather than on other grains of the top layer). Hence the lower layers provide for the top one a washboard potential, as shown in Fig. 1.

We place the grains of the top layer on the sites of a regular square lattice with row index $`t`$ and columns $`i`$. At any given time a grain $`G`$ may become active if at least one of its neighbors from the row above has been active at the previous time step. If $`\mathrm{\Delta }E(G)`$, the total energy transferred from these neighbors, exceeds the barrier $`E_b`$ of the washboard, $`G`$ becomes active, “rolls down” and collides with the grains of the next row. The energy it brings to these collisions is $`1+\mathrm{\Delta }E(G)`$, where 1 is our unit of energy, representing the potential energy due to the height difference between two consecutive rows (see Fig. 1). A fraction $`f`$ of its total energy is dissipated; the rest is divided stochastically among its three neighbors from the lower row.

The model is defined in terms of two variables; an activation variable $`S_i^t=0,1`$ and an energy $`E_i^t`$. The index $`t`$ denotes rows of our square lattice and time; at time $`t`$ we update the states of the grains that belong to row $`t`$. The model is controlled by two parameters: $`E_b`$, the barrier height, and $`f`$, the fraction of dissipated energy.

The dynamic rules of our model are as follows. For given activities $`S_i^t`$ and energies $`E_i^t`$ we first calculate the energy transferred to the grains of the next row $`t+1`$. To this end we generate for each active site three random numbers, $`z_i^t(\delta )`$ (with $`\delta =\pm 1,0`$) that add up to 1. The energy transferred to grain $`(t+1,i)`$, given by

$$\mathrm{\Delta }E_i^{t+1}=(1f)\underset{\delta =\pm 1,0}{}S_{i\delta }^tE_{i\delta }^tz_{i\delta }^t(\delta ),$$
(4)   
determines its activation:

$$S_i^{t+1}=\{\begin{array}{cc}1& \text{active \hspace{0.33em}\hspace{0.33em} if }\mathrm{\Delta }E_i^{t+1}>E_b,\\ 0& \text{inactive if }\mathrm{\Delta }E_i^{t+1}E_b.\end{array}$$
(5)   
Then the energies of the next row of grains are set:

$$E_i^{t+1}=S_i^{t+1}(1+\mathrm{\Delta }E_i^{t+1}).$$
(6)   
The three random numbers $`z_i^t(\delta )`$ represent the fraction of energy transferred from the grain at site $`(t,i)`$ to the one at $`(t+1,i+\delta )`$. We add up the energy contributions from these active sites; the fraction $`1f`$ is not dissipated; if the acquired energy $`\mathrm{\Delta }E_i^{t+1}`$ exceeds $`E_b`$, site $`(t+1,i)`$ becomes active, rolls over the barrier and brings to the collisions (at time $`t+2`$) the acquired energy calculated above and its excess potential energy (of value 1).

## 3 Qualitative discussion of the transition and connection to the experiment.

Let us vary $`E_b`$ at a fixed value of the dissipation. For small values of $`E_b`$ an active grain will activate the grains below it with high probability; avalanches will propagate downhill and also spread sideways. For a strongly localized initial activation we should, therefore, see triangular shaped activated regions. As $`E_b`$ increases, the rate of activation decreases and the opening angle $`\theta `$ of these triangles should decrease, until $`E_b`$ reaches a critical value $`E_b^c`$, beyond which initial activations die out in a finite number of time steps (or rows). These expectations are indeed borne out by simulations of the model: the dependence of $`E_b^c`$ on the dissipation $`f`$ is shown in Fig. 2.

The physics of the process is captured by a simple mean-field type approximation, in which all stochastic variables are replaced by their average values. Consider an edge separating an active region from an inactive one. At time $`t`$ sites to the left of $`i`$ and $`i`$ itself are wet, whereas $`i+1,i+2,\mathrm{}`$ are dry. Will site $`i`$ be wet or dry at the next time step? In our mean-field estimate of the answer, assuming that all wet sites at time $`t`$ have the same energy $`E^t`$, the energy delivered to site $`i`$ at time $`t+1`$ is $`\mathrm{\Delta }E_i^{t+1}=\frac{2}{3}(1f)(1+\mathrm{\Delta }E^t)`$, where we set in Eq. (4) all $`z(\delta )=1/3`$. At the critical point we expect all energies to just suffice to go over the barrier; hence set $`\mathrm{\Delta }E_i^{t+1}=\mathrm{\Delta }E^t=E_b^c`$. Solving the resulting equation yields a rough estimate of the transition line,

$$E_b^c=2(1f)/(1+2f),$$
(7)   
as shown in Fig. 2. It is easy to produce better mean-field type estimates of the transition and to compute the corresponding energy profile of the wet region .

To connect our model to the DD experiments note that the tilt angle $`\phi `$ tunes the ratio between the barrier height and the difference of potential energies between two rows. When the system is prepared at $`\phi _0`$, this ratio is precisely $`E_b^c`$. When one increases the tilt angle to $`\phi >\phi _0`$, $`E_b`$ (measured in units of the potential difference) decreases and we have $`E_b<E_b^c`$. As the tilt angle is now reduced, the size of $`E_b`$ increases, until it reaches its critical value precisely at $`\phi _0`$. Thus increasing $`E_b`$ in the model corresponds to lowering the tilt angle towards $`\phi _0`$ where the system is precisely at its boundary of dynamic stability.

Hence to reproduce the experiment we were looking for (a) fairly compact triangular regions of activation for $`E_b<E_b^c`$, and (b) an opening angle of these triangles which should go to zero as $`E_b`$ approaches $`E_b^c`$ from below.

We simulated the model defined in Eqs. (4)-(6) and found that it indeed reproduces these qualitative features of the experiment (see Fig. 3). The avalanches shown were produced for dissipation $`f=0.5`$, activating a single site at $`t=0`$, to which an initial energy of $`E_0=100`$ was assigned. As long as $`E_b`$ was not too close to $`E_b^c`$ the observed avalanches were compact, triangular and with fairly straight edges. The edges became rough only very close to $`E_b^c`$, such as the one shown on the right hand side of Fig. 3. The opening angle of the active regions $`\theta `$ decreased as $`E_b`$ increased towards $`E_b^c`$, as indicated in in the inset of Fig. 2. From these simulations we estimated $`E_b^c`$ and the exponent (see Eq. (3))

$$x=\nu _{}\nu _{}=0.98(5)1.$$
(8)   
The linear variation of $`\mathrm{tan}(\theta )`$ with $`\mathrm{\Delta }\phi `$ is in agreement with experimental measurements . Our findings have to be compared with the mean-field theory suggested in Ref.  which predicts a square root behavior.

## 4 Crossover to directed percolation.

The linear law (8) is consistent with the critical exponents of CDP 

$$\nu _{}=2,\nu _{}=1,\beta =0.$$
(9)   
These observations pose, however, a puzzle: since one believes that DP is the generic situation, one would expect to find non-compact active regions and DP exponents. In fact, according to the DP conjecture  any continuous spreading transition from a fluctuating active phase into a single frozen state should belong to the universality class of directed percolation (DP), provided that the model is defined by short range interactions without exceptional properties such as higher symmetries or quenched randomness. The present model has neither special symmetries nor randomness ; it has a fluctuating active phase and exhibits a transition, characterized by a positive one-component order parameter, into a single absorbing state. Hence the phase transition of our model should belong to the DP universality class.

In order to understand this apparent paradox we performed high-precision Monte-Carlo simulations for dissipation $`f=0.5`$ (see  for further details). We performed time-dependent simulations , i.e., we toppled a single grain in the center of the top row and measured the survival probability $`P(t)`$ and the number of active sites $`N(t)`$. At criticality, these quantities exhibit an asymptotic power law behavior

$$P(t)t^\delta ,N(t)t^\eta .$$
(10)   
In the case of CDP these exponents are given by  $`\delta =1/2`$ and $`\eta =0`$, whereas DP is characterized by the exponents  $`\delta =0.1595`$ and $`\eta =0.3137`$. Detecting deviations from power-law behavior in the long-time limit we estimated $`E_b^c=0.385997(5)`$.

Numerical results, obtained from simulations at $`E_b^c`$, are shown in Fig. 4. After a short transient the system enters an intermediate regime, which extends up to several hundred time steps. Here the active sites form a single cluster and we observe power-law behavior with CDP exponents (dotted lines in Fig. 4). This intermediate regime is followed by a long crossover from CDP to DP, extending over almost two decades up to more than $`10^4`$ time steps, after which the system enters an asymptotic DP regime (indicated by dashed lines in Fig. 4).

Compared with ordinary DP lattice models, this crossover regime is extremely long. We observed that by increasing $`f`$ the crossover time can be reduced by more than one decade. Hence, for an experimental verification of DP, systems with high dissipation are more appropriate. The present experiments correspond to about 3000 time steps (rows of beads); increasing this to about $`10^4`$ by using a longer inclined plane and smaller beads should yield DP behavior, provided that deviations of the experiment from the model do not increase with system size.

The crossover from CDP to DP is illustrated in Fig. 5. Two avalanches are plotted on different scales. The left one represents a typical avalanche within the first few thousand time steps. As can be seen, the cluster appears to be compact on a lateral scale up to 100 lattice sites. However, as can be seen in the right panel of Fig. 5, after very long time the cluster breaks up into several branches, displaying the typical patterns of critical DP clusters. Thus, before measuring critical exponents, this feature has to be tested experimentally. To this end the DD experiment should be performed repeatedly at the critical tilt $`\phi =\phi _0`$. In most cases the avalanches will be small and compact. However, sometimes large avalanches will be generated which reach the bottom of the plate. If these avalanches are non-compact, we expect DP-type asymptotic critical behavior. Only then is it worthwhile to optimize the experimental setup and to measure the critical exponents quantitatively.

AJD thanks the CICPB and the UNLP, Argentina and the Weizmann Institute for financial support. HH thanks the Weizmann Institute and the Einstein Center for hospitality and financial support. ED thanks the Germany-Israel Science Foundation for support, B. Derrida for some most helpful initial insights and A. Daerr for communicating his results at an early stage.