Anchored advected interfaces, Oslo model, and roughness at depinning

Here we introduce a new model, the DEP. This is the simplest model we found that breaks left-right symmetry, equivalent to advection. Depending on whether particle-hole symmetry is unbroken or broken, it belongs to the Edwards–Wilkinson or KPZ class.

In this model, each site $x = 1,2,\dots$ is either empty or occupied. Particles jump with rate 1 to an empty neighbour as in the simple exclusion process, but an additional directed transition is introduced, which allows particles to jump over a particle to its right (a jump of length 2) with rate 1, and holes over a hole to the right (a jump of length 2) with rate r. The transitions are

• 1.  
  10 → 01 with rate 1,
• 2.  
  01 → 10 with rate 1,
• 3.  
  110 → 011 with rate 1,
• 4.  
  001 → 100 with rate r.

A rate r = 1 preserves the particle-hole symmetry, while breaking the directional (left-right) symmetry. We refer to this model as the symmetric DEP (sDEP). The case r ≠ 1 is referred to as the asymmetric DEP (aDEP), which breaks both directional and particle-hole symmetry.

When using periodic boundary conditions, the product measure is stationary. This can be derived from the (non-detailed) balance equation. Then the total particle current $j(\rho)$ is given by

$$\begin{equation} j(\rho) = 2 \left[\rho^2(1-\rho) - r (1-\rho)^2 \rho \right]= 2\rho(1-\rho)\left[ \rho - r(1-\rho)\right]. \end{equation} \tag{ 3 }$$

Therefore, $\rho_{\mathrm{c}} = \frac{r}{1+r}$ is the critical density, where the overall current vanishes.

In contrast, we consider a large system with closed boundaries, where no particle current is allowed. This leads to self-organized criticality, where particle (or hole) excess will be pushed away to the right. Taking $L\to \infty$ results in a stationary density $\rho = \rho_{\mathrm{c}} = \frac{r}{1+r}$ in the domain of observation.

Toom's model is the subject of many theoretical and numerical studies [24, 29–33] ⁴ . As in the DEP, each site is either empty or occupied, but while in the DEP the jump range is either 1 or 2, in Toom's model it is unbounded. More precisely, one chooses a spin i, and exchanges it with the next spin to its right which has the opposite sign. The rates, for any $k\in \mathbb{N}$, are

• (1)k  
  $\underbrace{1\dots 1}_{k}0 \to 0\underbrace{1\dots 1}_k$ with rate 1,
• (2)k  
  $\underbrace{0\dots 0}_{k}1 \to 1\underbrace{0\dots 0}_k$ with rate r.

When r = 1, as in the DEP, the particle-hole symmetry is preserved, and we refer to the model as the symmetric Toom model (sToom). When r ≠ 1 this symmetry is broken, and we refer to the model as the asymmetric Toom model (aToom).

In Toom's model with periodic boundary conditions each stationary state is non-interacting, i.e. occupations at different sites are independent. The total particle current is

$$\begin{align} j(\rho)&=\rho (1-\rho) + 2 \rho^2 (1 - \rho) + 3 \rho^3(1-\rho)+ \dots - r \left[(1-\rho)\rho + 2(1-\rho)^2 \rho + \dots \right] \nonumber\\ &= \frac{\rho}{1-\rho}-\frac{r(1-\rho)}{\rho}. \end{align} \tag{ 4 }$$

This current vanishes at the critical density $\rho_{\mathrm{c}} = \frac{\sqrt r}{1+\sqrt{r}}$. As thoroughly discussed in [24], the system with closed boundaries and for $L\to \infty$ reaches a stationary state with this density at its left end.

As suggested in [24, 32], these models can be described using a continuum theory. We define the height function $h(x,t)$ as the number of particles between the left boundary (first spin at site 1) and x, minus its expectation. In the continuum limit we expect it to solve a stochastic differential equation of the type

$$\begin{align} \partial_t h(x,t) = c + D \nabla^2 h(x,t) -\mu \nabla h(x,t) + \lambda_2 \big[\nabla h(x,t)\big]^2 + \lambda_3 \big[\nabla h(x,t)\big]^3 + \dots + \sigma \xi(x,t).\nonumber\\ \end{align} \tag{ 5 }$$

Discarding higher-order terms, we are left with two possible limits: for sDEP and sToom, the particle-hole symmetry forces the equation to stay invariant under the transformation $h \mapsto -h$, hence the first and fourth terms are not allowed ($c = \lambda_2 = 0$), and we are left with the anchored advected Edwards–Wilkinson equation

$$\begin{equation} \begin{split} \partial_t h(x,t) =&\; D \nabla^2 h(x,t) - \mu \nabla h(x,t) + \sigma \xi(x,t),\\ h(t,0) =&\; 0 , \quad \left\langle\xi(x,t) \xi(x^{^{\prime}},t^{^{\prime}})\right\rangle = \delta(x-x^{^{\prime}})\delta(t-t^{^{\prime}}). \end{split} \end{equation} \tag{ aaEW }$$

In aDEP and aToom the constant and quadratic terms are allowed, leading to the anchored advected KPZ equation

$$\begin{equation} \begin{split} \partial_t h(x,t) =&\; D \nabla^2 h(x,t) -\mu \nabla h(x,t) + \lambda_2 \big[\nabla h(x,t)\big]^2 + c + \sigma \xi(x,t),\\ h(t,0) =&\; 0, \quad \left\langle\xi(x,t) \xi(x^{^{\prime}},t^{^{\prime}})\right\rangle = \delta(x-x^{^{\prime}})\delta(t-t^{^{\prime}}). \end{split} \end{equation} \tag{ aaKPZ }$$

Since $\left\langle h(x,t)\right\rangle = 0$ by construction, the first two terms on the rhs vanish on average, and hence

$$\begin{equation} c = -\lambda_2 \langle [\nabla h]^2 \rangle. \end{equation} \tag{ 6 }$$

An additional term with $\lambda_3\neq 0$ may lead to logarithmic corrections [34].

A direct derivation of this limit for Toom's model appears in [24, 32]. We describe here the derivation for the DEP.

Let us consider the DEP at the critical density $\rho_c = \frac{r}{1+r}$. Define $\eta(x,t)$ to be 1 if there is a particle at x at time t and 0 otherwise. Then the height function h is defined as $h(x,t): = \sum_{y = 1}^x [\eta(y,t) - \rho]$.

During a time period $\mathrm{d}t$, the value of h at x changes when particles jump from the left of x to its right or vice versa:

• 1.  
  $h(x,t+\mathrm{d}t) = h(x,t) - 1$ with probability $p_1 = \eta(x,t)[1-\eta(t,x+1)]\mathrm{d}t$,
• 2.  
  $h(x,t+\mathrm{d}t) = h(x,t) + 1$ w.p. $p_2 = [1-\eta(x,t)]\eta(t,x+1)\mathrm{d}t$,
• 3.  
  $h(x,t+\mathrm{d}t) = h(x,t) - 1$ w.p. $p_3 = \eta(x,t)\eta(x+1,t)[1-\eta(x+2,t)]\mathrm{d}t$,
• 4.  
  $h(x,t+\mathrm{d}t) = h(x,t) - 1$ w.p. $p_4 = \eta(x-1,t)\eta(x,t)[1-\eta(x+1,t)]\mathrm{d}t$,
• 5.  
  $h(x,t+\mathrm{d}t) = h(x,t) + 1$ w.p. $p_5 = r[1-\eta(x,t)][1-\eta(x+1,t)]\eta(x+2,t) \mathrm{d}t$,
• 6.  
  $h(x,t+\mathrm{d}t) = h(x,t) + 1$ w.p. $p_6 = r[1-\eta(x-1,t)][1-\eta(x,t)]\eta(x+1,t) \mathrm{d}t$.

First, we calculate

$$\begin{equation*} \langle h(x,t+\mathrm{d}t) - h(x,t) \rangle = (-p_1 + p_2 -p_3- p_4+p_5+p_6 ) \mathrm{d}t. \end{equation*}$$

We separate $\eta(x,t)$ into its average ρ plus fluctuations, $\eta(x,t) = \rho + \nabla h(x,t)$. This yields

$$\begin{align} \langle h(x,t+\mathrm{d}t) - h(x,t) \rangle &= \big[ D \nabla^2 h(x,t) - \mu \nabla h(x,t) + \lambda_2 \big(\nabla h(x,t)\big)^2 + c\big] \mathrm{d}t, \end{align} \tag{ 7 }$$

$$\begin{align} \qquad\qquad~ D &= 1+3r-6r\rho +3r\rho^2 + 3\rho^2, \end{align} \tag{ 8 }$$

$$\begin{align} \qquad\qquad~\; \mu &= -6 r \rho ^2+8 r \rho -2 r -6 \rho ^2+4 \rho \end{align} \tag{ 9 }$$

$$\begin{align} \qquad\qquad\;\, \lambda_2 &= 6 r\rho - 4r + 6\rho - 2. \end{align} \tag{ 10 }$$

The noise term is obtained from

$$\begin{equation} \sigma^2 \mathrm{d}t= \left\langle \frac{1}{k} \left[ \sum_{x=1}^k h(x,t+\mathrm{d}t)-h(x,t) \right]^2 \right\rangle. \end{equation} \tag{ 11 }$$

The reason we take a spatial average is to account for correlations: since each of the transitions (1) and (2) contribute 1 to the sum and transition (3) which is equivalent to (4), and transition (5) which is equivalent to (6), contribute 2, we obtain

$$\begin{equation} \left\langle \frac{1}{k} \left( \sum_{x=1}^k h(x,t+\mathrm{d}t)-h(x,t) \right)^2 \right\rangle = p_1+p_2+4p_3+4p_5. \end{equation} \tag{ 12 }$$

Note that $p_4,p_6$ are spatial shifts of $p_3,p_5$ hence we must count them only once. This yields

$$\begin{equation} \sigma^2 = 2 \rho(1-\rho) (1 +2 r +2 \rho-2 r \rho). \end{equation} \tag{ 13 }$$

Combining these two equations, we obtain for $\rho = \rho_c = \frac{r}{1+r}$

$$\begin{equation} D= \frac{4 r +1}{r +1}, \quad \mu = \frac{2 r }{r +1}, \quad \lambda_2 = 2(r-1), \quad \sigma^2= \frac{2 r (5 r +1)}{(r +1)^3}. \end{equation} \tag{ 14 }$$

This is model (aaEW) when r = 1 and (aaKPZ) when r ≠ 1. For the symmetric case this gives

$$\begin{equation} D=\frac52, \quad \mu = 1, \quad \lambda_2=0, \quad \sigma^2=\frac32. \end{equation} \tag{ 15 }$$

As discussed for Toom's model and the DEP, models in the advected Edwards–Wilkinson universality class have a stationary measure which depends on the boundary condition. A similar behaviour is observed in the continuum limit. If h is a solution of the diffusion equation (aaEW) with periodic boundary conditions, then $\tilde{h}(x,t) = h(x-\mu t,t)$ solves the non-advected Edwards–Wilkinson equation,

$$\begin{align} \partial_t \tilde{h}(x,t) &= \partial_t h(x-\mu t,t) - \mu \nabla h(x -\mu t,t) = D \nabla^2 h(x-\mu t,t) + \sigma \xi(x-\mu t,t) \nonumber\\ &= D\nabla^2 \tilde{h}(x,t)+\sigma {\xi}(x,t), \quad \left\langle\xi(x,t) \xi(x^{^{\prime}},t^{^{\prime}})\right\rangle = \delta(x-x^{^{\prime}})\delta(t-t^{^{\prime}}). \end{align} \tag{ 16 }$$

(We used that $\xi(x,t)$ and $\xi(x-\mu t,t)$ have the same correlations.) In particular, the stationary state for a periodic system of size L is a Brownian motion with periodic boundary conditions (i.e. a Brownian bridge). For $L\to \infty$, this change of variables is valid on the entire line $x\in\mathbb{R}$. However, if we consider the half line $x \in [0,\infty)$ with Dirichlet boundary condition $h(0) = 0$ as in equation (aaEW), we cannot define $\tilde{h}$ as above, and the stationary state is not a Brownian. This choice of boundary, corresponding to an anchored interface, is the subject of this article.

When considering portions of the interface far away from the origin, the effect of the boundary becomes negligible, and in particular the stationary state on an interval $[x,x+\Delta x]$ looks like a Brownian motion if $\Delta x \ll x$. This is shown rigorously for Toom's model in [33, 35], and we expect a similar behaviour in all models of this universality class.

As mentioned above, unlike the bulk behaviour, the behaviour of the model near the edge is drastically different from the standard Edwards–Wilkinson universality class [24, 28, 32]. The combination of the Dirichlet boundary and the advection term transfers temporal correlations into special ones. In the following, we derive analytic results for the edge behaviour in the anchored case.

We solve the advected diffusion equation with Dirichlet boundary conditions, see also [28] ⁵ . By simple rescaling we write equation (aaEW) in the form:

$$\begin{equation} \begin{split} &\partial_t h(x,t) = D\nabla^2 h(x,t) - \mu \nabla h(x,t) + \xi (x,t),\qquad h(x=0,t) = 0.\\ & \left\langle \xi (x,t) \xi(x^{^{\prime}},t^{^{\prime}})\right\rangle = \delta(x-x^{^{\prime}}) \delta(t-t^{^{\prime}}). \end{split} \end{equation} \tag{ 17 }$$

Equation (17) for x > 0 is solved by

$$\begin{equation} h(x,t) = \int_{y\gt 0} \int_{t^{^{\prime}}\lt t} P^{\mu}(x,t|y,t^{^{\prime}}) \xi(y,t^{^{\prime}}). \end{equation} \tag{ 18 }$$

Here

$$\begin{align} P_0(x,t) &= \frac{e^{-\frac{x^2}{4D t}}}{\sqrt{4 D\pi t}}, \nonumber \\ P^\mu_0(x,t) &= P_0(x,t) \,\mathrm{e}^{\frac {\mu x}{2D}-\frac{\mu^2 t}{4D}},\nonumber\\ P^\mu(x,t|y,t^{^{\prime}}) &= \big[ P_0 (|x-y|,t-t^{^{\prime}})-P_0 (x+y,t-t^{^{\prime}})\big] \,\mathrm{e}^{\frac {\mu (x-y)}{2D}-\frac{\mu^2 (t-t^{^{\prime}})}{4D}}. \end{align} \tag{ 19 }$$

In this section we calculate $\left\langle h(x,t)^2\right\rangle$ analytically. The reader only interested in the large-scale behaviour may skip to the next section.

The roughness exponent ζ is given by the behaviour of $\left\langle h(x,t)^2\right\rangle$ in the steady state,

$$\begin{align} \left\langle h(x,t)^2\right\rangle &= \int_{y_1\gt 0} \int_{t_1\lt t} P^{\mu}(x,t|y_1,t_1) \int_{y_2\gt 0} \int_{t_2\lt t} P^{\mu}(x,t|y_2,t_2) \left\langle\xi(y_1,t_1) \xi(y_2,t_2)\right\rangle \nonumber\\ &= \int_{y\gt 0} \int_{t^{^{\prime}}\lt t} P^{\mu}(x,t|y,t^{^{\prime}})^2 . \end{align} \tag{ 20 }$$

This is independent of t, and we will drop t from now on. To proceed, we use the Fourier-transform of the propagator (19),

$$\begin{align} P^{\mu}(x,t|y,0) &= \frac{\mathrm{e}^{-\frac{\mu ^2 t}{4D}+\frac{\mu (x-y)}{2D} }}{ \sqrt{4 D \pi t}} \left[\mathrm{e}^{- \frac{(x-y)^2}{4 Dt}} -\mathrm{e}^{-\frac{(x+y)^2}{4 Dt}} \right] \nonumber\\ &=\int \frac{\mathrm{d k}}{2\pi} \left[\mathrm{e}^{i k (x-y)}-\mathrm{e}^{i k (x+y)}\right] \mathrm{e}^{-k^2 D t-\frac{\mu ^2 t}{4D}+\frac{\mu }{2D} (x-y)} . \end{align} \tag{ 21 }$$

From now on, we set $\mu = D = 1$. We can recover the full µ, D and σ-dependence by remarking that

$$\begin{equation} \left\langle h(x)^2\right\rangle_\mu = \frac{\sigma^2}\mu \left\langle h\left(\frac{\mu x}{D},\frac{\mu^2 t}D\right)^2\right\rangle_{\mu=D=\sigma=1}. \end{equation} \tag{ 22 }$$

Then

$$\begin{equation} \int_{y\gt 0}\int_{t\gt 0} P^{\mu=1}(x,t|y,0)^2 = \int_{k,p}\frac{-16 \, k p \mathrm{e}^{x [1 +i (k+p)]}}{[2(k^2+p^2 )+1] [(k-p)^2+1 ] [(k+p)^2+1 ] }. \end{equation} \tag{ 23 }$$

Integrating equation (23) over k yields

$$\begin{equation} \left\langle h(x)^2\right\rangle_{\mu=1} = \int_{-\infty}^{\infty}\frac {\mathrm{d} p}{2\pi i} \frac{16 p \mathrm{e}^{-\sqrt{p^2+\frac{1}{2}} x+i p x+x}}{(4 p^2+1 )^2} . \end{equation} \tag{ 24 }$$

Equation (24) can be simplified, using for the denominator

$$\begin{equation} \int_0^\infty \mathrm{d} t\, t \,\mathrm{e}^{-a t} = \frac1{a^2}, \end{equation} \tag{ 25 }$$

and for $\mathrm{e}^{-\sqrt{p^2+\frac{1}{2}} x}$

$$\begin{equation} \int_{s\gt 0} \frac{\mathrm{e}^{-a s-\frac{1}{4 s}}}{2 \sqrt{\pi } s^{3/2}} =\mathrm{e}^{-\sqrt a}. \end{equation} \tag{ 26 }$$

Then the integral over p can be done, leading to

$$\begin{equation} \left\langle h(x)^2\right\rangle_{\mu=1} = \int_{s\gt 0}\int_{t\gt 0} \frac{2 t x \,\mathrm{e}^{-\frac{x^2}{4 s x^2+16 t}-\frac{s x^2}{2}-\frac{1}{4 s}-t+x}}{\pi s^{3/2} (s x^2+4 t )^{3/2}}. \end{equation} \tag{ 27 }$$

Setting $s \to 4st/x^2$, and integrating over t yields

$$\begin{equation} \left\langle h(x)^2\right\rangle_{\mu=1} = \int_{s\gt 0}\frac{\mathrm{e}^x x^2 K_0\left(\frac{(2 s+1) x}{2 \sqrt{s} \sqrt{s+1}}\right)}{4 \pi s^{3/2} (s+1)^{3/2}} . \end{equation} \tag{ 28 }$$

K₀ is the Bessel K₀-function. This can further be simplified, by first setting $y = \sqrt{s/({1+s})}$, and second $2u = y+y^{-1}$. The final result is

$$\begin{equation} \left\langle h(x)^2\right\rangle_{\mu=1} =\frac {x^2 \mathrm{e}^{x}}{\pi} \int_1^{\infty} K_0(x u)\, {\mathrm{d} u} = \frac {x \mathrm{e}^{x}}{\pi} \int_x^{\infty} K_0(u)\, {\mathrm{d} u} . \end{equation} \tag{ 29 }$$

The integral is known analytically,

$$\begin{equation} \left\langle h(x)^2\right\rangle_{\mu=1} = \frac{x \,\mathrm{e}^x }{2} \Big[ 1-x \left(\pmb{L}_{-1}(x) K_0(x)+\pmb{L}_0(x) K_1(x) \right) \Big]. \end{equation} \tag{ 30 }$$

The function $\pmb{L}_n(x)$ is the modified Struve-function. For µ ≠ 1 this reads according to equation (22)

$$\begin{equation} {\left\langle h(x)^2\right\rangle =\frac {x \,\mathrm{e}^{\mu x}}{\pi} \int_{\mu x}^{\infty} K_0(u)\, {\mathrm{d} u} } . \end{equation} \tag{ 31 }$$

One can approximate the Bessel function for large argument as

$$\begin{equation} K_0(x) = \mathrm{e}^{-x}\left[ \sqrt{\frac{\pi }{2 x}}+\mathcal{O}(x^{-\frac 32}) \right]. \end{equation} \tag{ 32 }$$

Then the integration can be done analytically,

$$\begin{equation} \left\langle h(x)^2\right\rangle \simeq \frac{e^{\mu x} x\, \text{erfc}\left(\sqrt{\mu x}\right)}{\sqrt {2}} = \frac{\sqrt{x}}{\sqrt{2 \pi\mu }} + \mathcal{O}(x^{-\frac12}). \end{equation} \tag{ 33 }$$

This gives a roughness exponent

$$\begin{equation} \zeta= \frac14. \end{equation} \tag{ 34 }$$

From figure 1 we see that, for large x, we can approximate equation (18) as

$$\begin{align} h(x,t=0) &= \int_0^\infty \mathrm{d} y \int_{-\infty}^{0}\mathrm{d} t^{^{\prime}} P^{\mu}(x,t|y,t^{^{\prime}}) \xi(y,t^{^{\prime}}) \nonumber\\ &\approx \int_{-\infty}^\infty \mathrm{d} y \int^0_{-\frac x\mu } \mathrm{d} t^{^{\prime}}\,P_0^\mu(x-y,-t^{^{\prime}}) \xi(y,t^{^{\prime}})\nonumber\\ &= \int_{-\infty}^\infty \mathrm{d} y \int_0^{\frac x\mu } \mathrm{d} t^{^{\prime}}\,P_0^\mu(x-y,t^{^{\prime}}) \xi(y,t^{^{\prime}}). \end{align} \tag{ 35 }$$

Then

$$\begin{align} \left\langle h(x)^2\right\rangle &\approx \int_{-\infty}^\infty \mathrm{d} y \int_0^{\frac x\mu } \mathrm{d} t^{^{\prime}}\,P_0^\mu(x-y,t^{^{\prime}})^2 \nonumber\\ &= \frac{1}{\sqrt{8\pi}} \int_0^{\frac x\mu } \frac{\mathrm{d} t^{^{\prime}}}{\sqrt{t^{^{\prime}}}} = \sqrt{\frac{x}{2\pi \mu }}. \end{align} \tag{ 36 }$$

This reproduces the result of equation (33).

Zoom In Zoom Out Reset image size

We wish to calculate (with $\mu = D = 1$)

$$\begin{equation} \left\langle h(x,t) h(x,0) \right\rangle = \int_{y\gt 0} \int_{t^{^{\prime}}\lt 0} P_{\mu}(x,0|y,t^{^{\prime}}) P_{\mu}(x,t|y,t^{^{\prime}}) . \end{equation} \tag{ 37 }$$

With the help of equation (35), this can be approximated as

$$\begin{align} \left\langle h(x,t) h(x,0) \right\rangle \approx \int_y \int_{0\leqslant t^{^{\prime}}\leqslant x} P_{0}^\mu(y,t^{^{\prime}}) P_{0}^\mu(y+t,t+t^{^{\prime}}). \end{align} \tag{ 38 }$$

Therefore

$$\begin{align} \!\left\langle h(x,t) h(x,0) \right\rangle &\simeq \int_{-\infty}^\infty \frac{\mathrm{d} k}{2\pi}\int_{-\infty}^\infty \frac{\mathrm{d} p}{2\pi} \int_y \int_{0\leqslant t^{^{\prime}}\leqslant x} \mathrm{e}^{ik(x-y+t^{^{\prime}}) -k^2t^{^{\prime}}} \mathrm{e}^{i p (x-y+t+t^{^{\prime}}) -p^2 (t+t^{^{\prime}})}\nonumber\\[4pt] \quad& = \int_{0\leqslant t^{^{\prime}}\leqslant x} \int_{-\infty}^\infty \frac{\mathrm{d} k}{2\pi} \mathrm{e}^{-k^2 (t+2t^{^{\prime}}) -i k t } =\int_{0\leqslant t^{^{\prime}}\leqslant x} \frac{\mathrm{e}^{-\frac{t^2}{4(t+2t^{^{\prime}})}}}{2 \sqrt{\pi} \sqrt{t+2 t^{^{\prime}}}} \nonumber\\[4pt] \quad&=\frac1 4 \sqrt{\frac t \pi}\int_{0}^{\frac {2 x }t}\mathrm{d} y\, \frac{\mathrm{e}^{-\frac{t}{4(1+y)}}}{\sqrt{1+y}} = \frac t{8 \sqrt{\pi} }\left[ \Gamma\left(-\frac12,\frac{t^2}{4(t+2x)} \right) -\Gamma\left(-\frac12,\frac{t}{4} \right)\right] . \end{align} \tag{ 39 }$$

The Γ-function decays as a power-law corrected Gaussian for large second argument. As we are interested in both large x and large t, the second term vanishes; the first one implies that the proper scaling variable is $t/\sqrt x$, reducing the denominator of the first term from $t+2x\to 2x$. As a result,

$$\begin{align} \left\langle h(x,t) h(x,0) \right\rangle &\simeq \sqrt{\frac{{x}}{{2\pi}}} \, f\!\left(\frac t{\sqrt x}\right), \quad f(0)=1, \quad f(y) \stackrel{y\to \infty}{\longrightarrow} 0, \nonumber\\[4pt] f (y) &= \frac{y }{4 \sqrt{2}} \,\Gamma \left(-\frac{1}{2},\frac{y^2}{8}\right). \end{align} \tag{ 40 }$$

We tested this against a numerical integration of the exact expression. This worked, but only for large x and t; for small times our approximation converges from the wrong side. f(y) has series expansions for small and large y

$$\begin{align} f(y) &= 1-\frac{1}{2} \sqrt{\frac{\pi }{2}} y+\frac{y^2}{8}-\frac{y^4}{384}+\frac{y ^6}{15\,360} + \cdots \nonumber\\[4pt] f(y) &= \mathrm{e}^{-\frac{y^2}8}\left[\frac{4}{y^2}-\frac{48}{y^4}+\frac{960}{y^6}+ \cdots \right] . \end{align} \tag{ 41 }$$

Putting back the dependence on µ, D and σ yields with equation (22)

$$\begin{equation} \left\langle h(x,t) h(x,0) \right\rangle \simeq \sigma^2 \sqrt{\frac{x}{2\pi \mu D}}\, f\!\left(\frac{\mu^2 t}{\sqrt{D x \mu}}\right). \end{equation} \tag{ 42 }$$

Our calculations above teach us that

$$\begin{equation} z= \left \{ \begin{array}{cl} 2& \textrm{for bulk observables in the comoving (advected) frame} \\[4pt] 1& \textrm{for bulk observables in the fixed frame} \\[4pt] \frac12 & \textrm{for } \left\langle h(x,t) h(x,0) \right\rangle \end{array} \right. \end{equation} \tag{ 43 }$$

In particular, we stress that at a time scale $t\sim L$ we reach stationarity on the interval $[0,L]$ (see also [33, 35]).

In the symmetric case, we expect the (centred) height function h in both Toom's model and the DEP to converge to the solution of the aaEW equation (aaEW) (up to possible logarithmic corrections). In order to test this, we compare the correlation $\langle h(x,t)h(x,0) \rangle$ to equation (42).

In figure 2 we show this comparison for Toom's model. It is worth noting, as already mentioned in [24], that there are no finite-size effects in this simulation. We see that the analytic prediction (42) is well verified. (We did not check the dependence on parameters µ, D, and σ). For more thorough numerical studies of this model we refer the reader to [24, 32, 36].)

Zoom In Zoom Out Reset image size

In figure 3 we show the same plot for the sDEP. There are noticeable finite-size effects, but for large systems the rescaled correlator $\left\langle h(x,t) h(0,0\right\rangle x^{-1/2}$ plotted against $t/x^2$ seemingly converges. Remarkably, we obtain not only the correct exponent $\zeta = \frac{1}{4}$, but also the correct coefficient: after rescaling equation (17) to account for the coefficients D and σ in equation (aaEW), and plugging in $D,\mu$ and σ found in section 2.4, we obtain

$$\begin{align*} \langle h(x)^2 \rangle \approx \frac{\sigma^2}{\sqrt {2\pi \mu D}}\sqrt x = \sqrt{\frac{9}{20 \pi}} \sqrt x \approx 0.378 \sqrt x. \end{align*}$$

This relation, including the coefficient of 0.378 is verified numerically, as can be seen on figure 3. On the other hand, the rescaling factor for the argument of f is off by a factor of about 1.66,

$$\begin{equation} \left\langle h(x,t) h(0,0)\right\rangle x^{-1/2} \approx \sqrt{\frac{9}{20 \pi}} f \left( {1.66} \sqrt {\frac25}\frac t{\sqrt{x}} \right). \end{equation} \tag{ 44 }$$

This factor, which seemingly affects time scales only, may be due to the presence of sub-sub-leading corrections $\sim (\nabla h)^3$ in the equation of motion. The presence of this type of corrections has been discussed in [34, 36]. As a marginal perturbation $(\nabla h)^3$ gives logarithmic corrections, which may be of the form $\ln L$ or $\ln t$. Our simulation suggests that the combination $\frac{\sigma^2}{\sqrt{\mu D}}$ does not renormalize, while each coefficient by itself may renormalize. In particular the combination $\mu^3/D$ renormalizes by a factor of 2.8 at $L = 2^{13}$. We show in section 4.2 that $\sigma^2/D$ does not renormalize in the bulk. In summary, static observables seem not to renormalize, while temporal correlations do.

Zoom In Zoom Out Reset image size

The code used to simulate the DEP is available with this paper.

The exponent $z = \frac{1}{2}$ in the aaEW universality class can also be derived directly in the discrete model, without passing to the continuum limit.

Consider the height function h(L) of the sDEP at fixed position L. This function only changes when a particle or a hole at L or L + 1 jumps to the right of L. Since the density is approximately $1/2$, at any such jump h(L) increases by 1 with probability close to $1/2$ and decreases by 1 with probability close to $1/2$. However, due to density fluctuations, this probability is not exactly $1/2$. When the density is slightly above $1/2$ there are more particles, and h is more likely to decrease. Conversely, when the density fluctuates below $1/2$, it is more likely to increase. The exact correction is complicated, depending on non-trivial correlations, but to first order we may assume that it induces a drift $-\alpha h(L)$ with α > 0, possibly L-dependent. We conclude that h(L) behaves as a random walk with a drift term $-\alpha h(L)$, equivalent to diffusion in a confining potential $\frac{\alpha}{2} h(L)^2$. For this process, we know that h(L) fluctuates on a scale $\alpha^{-1/2}$ and that its relaxation time is of order $1/\alpha$. Finally, since $\zeta = \frac{1}{4}$, the fluctuations of h are of order $L^{1/4}$. Comparing this scaling with $\alpha^{-1/2}$ suggests that α scales as $L^{-1/2}$, and hence $z = \frac{1}{2}$.

The aaKPZ universality class can be analysed in a similar manner. As in the symmetric case, the stationary state in the periodic or infinite system is that of Brownian motion, and the bulk behaviour is determined by the KPZ equation in the comoving frame. The edge behaviour in the aaKPZ universality class can then be studied using methods close to those described above [24–27], yielding a roughness exponent $\zeta = \frac{1}{3}$, and dynamical exponents $z = \frac{3}{2}$ for bulk observables in the comoving frame, z = 1 for bulk observables in the fixed frame, and $z = \frac{2}{3}$ for $\langle h(x,t)h(x,0)\rangle$.

Consider the equation of motion (17), with a subleading KPZ-term, and a subsubleading term $\sim (\nabla h)^3$, ⁶

$$\begin{equation} \eta \partial_t h(x,t) = D\nabla^2 h(x,t) - \mu \nabla h(x,t) + \lambda_2 [\nabla h(x,t)]^2 + \lambda_3 [\nabla h(x,t)]^3 +\cdots+ \xi (x,t). \end{equation} \tag{ 45 }$$

In the absence of boundaries, the measure $\mathcal{P}_t[h]$ satisfies the Fokker–Planck equation

$$\begin{align} \eta {\partial _t}{{\cal P}_t}[h] = {\sigma ^2}\int_x {\frac{{{\delta ^2}}}{{\delta h{{(x)}^2}}}} {{\cal P}_t}[h] - \int_x {\frac{\delta }{{\delta h(x)}}} \left\{ {\left( {D{\nabla ^2}h(x) - \mu \nabla h(x,t) + \sum\limits_{n2} {{\lambda _n}} {{\left[ {\nabla h(x)} \right]}^n}} \right){{\cal P}_t}[h]} \right\}{\text{ }}. \end{align} \tag{ 46 }$$

Similar to what happens for KPZ, (see e.g. [2], section 7.12), for $\mu = \lambda_i = 0$, a steady-state solution $\partial_t \mathcal{P}^{\mathrm{ss}}[h] = 0$ can be found by asking that

$$\begin{equation} \sigma^2 \frac{\delta}{\delta h(x)} \mathcal{P}^{\mathrm{ss}}[h] = D \nabla^2 h(x) \mathcal{P}^{\mathrm{ss}}[h] . \end{equation} \tag{ 47 }$$

This is solved by

$$\begin{equation} \mathcal{P}^{\mathrm{ss}}[h] = \mathcal{N} \exp\left({-\frac{\sigma^2}{2D} \int_x [\nabla h(x)]^2} \right). \end{equation} \tag{ 48 }$$

What is the effect of the additional terms? Inserting the steady state into equation (46) yields

$$\begin{align} \eta\partial_t \mathcal{P}^{\mathrm{ss}}[h] &= - \int_x \frac{\delta}{\delta h(x)} \left\{\left( - \mu \nabla h(x,t)+ \sum_{n\geqslant 2} \lambda _n \left[ \nabla h(x) \right]^n \right) \mathcal{P}^{\mathrm{ss}}[h] \right\} \nonumber\\ &= \frac{\sigma^2}{D } \int_x \nabla^2 h(x) \left( - \mu \nabla h(x,t)+ \sum_{n \geqslant 2}\lambda _{n} \left[ \nabla h(x) \right]^n\right) \mathcal{P}^{\mathrm{ss}}[h] \nonumber\\ &= \frac{\sigma^2}{D } \int_x \nabla \left( - \frac\mu2 [\nabla h(x,t)]^2+ \sum_{n \geqslant 2} \frac{\lambda _n}{n+1} \left[ \nabla h(x) \right]^n \right) \mathcal{P}^{\mathrm{ss}}[h] =0. \end{align} \tag{ 49 }$$

Note that by going from the first to the second line, we have dropped the derivative of the terms inside the big round brackets, since they are a total derivative. The last line vanishes, again due to the fact that this is a total derivative.

This calculation shows that as long as we consider a system without boundaries, the steady state is independent of µ and λ_n . When we consider the system from an RG perspective, this means that $\sigma^2/D$ is not renormalized. There is nothing, however, to protect $\eta/D$. In the presence of a boundary, we do not expect bulk properties to renormalize differently, thus the effective η, D, µ, and λ_i to be unchanged if a boundary is introduced.

Let us mention, that from an RG perspective there may appear additional renormalizations on the boundary, and since $\left\langle h(x,t) h(x,0)\right\rangle$ depends on the distance x to the boundary, they could appear there. Since the static 2-point function $\left\langle h(x,t)^2 \right\rangle$ shows no corrections, and temporal-derivative terms are marginal in the bulk, we do not expect this to be the case. It would be interesting to render this intuitive argument more rigorously.

The Oslo model describes the evolution of a sand pile, given by its height h (number of grains) at each horizontal position (see figure 4). In a real sandpile, whether a grain at site i slides downhill depends on the local slope z(i), the friction with its neighbours of contact, and its orientation. The Oslo model is a simple one-dimensional model for this phenomenon, which depends only on the local slope z(i), and a random variable $z_{\mathrm{c}}(i)$. It was introduced in [22, 23], and is defined as follows: Consider the height function h(i) in the left of figure 4. To each height profile h(i) associate a stress field (slope) z(i) defined by

$$\begin{equation} z(i) := h(i)-h(i+1). \end{equation} \tag{ 50 }$$

Zoom In Zoom Out Reset image size

In addition, at each position there is a threshold $z_{\mathrm{c}}(i)$. A toppling is invoked if $z(i)\gt z_{\mathrm{c}}(i)$, i > 1. The toppling rules are

$$\begin{align} z(i) \to z(i)-2 , \quad z(i\pm 1)\to z(i\pm 1)+1 . \end{align} \tag{ 51 }$$

They can be interpreted as moving a grain from the top of the column at site i to the top of the column at site i + 1,

$$\begin{equation} h(i)\to h(i)-1,\quad h(i+1)\to h(i+1)+1. \end{equation} \tag{ 52 }$$

After such a move, the threshold $z_{\mathrm{c}}(i)$ for site i is updated,

$$\begin{align} z_{\mathrm{c}}(i) \to \textrm{new random number} . \end{align} \tag{ 53 }$$

In its original version, the random number is 1 or 2 with probability $1/2$. To obtain figure 4 we used a random number drawn uniformly from the interval $[0,2]$. This reduces the critical slope.

In this article we consider a sysmtem of size L with absorbing (Dirichlet) boundary conditions on the left boundary (site 1) and free (Neumann) boundary conditions on the right boundary (site L). That is, in the h variable, a toppling at 1 happens when $z(1)\gt z_{\mathrm{c}}(1)$, and it moves a grain from site 1 to site 2. On the right boundary, a toppling at L happens when $h(L)\gt z_{\mathrm{c}}(L)$, and it removes one grain from h. We can formulate this in the language of the z variables: when a toppling at 1 occurs, z(1) decreases by 2 and z(2) increases by 1 (i.e. the particle that should have gone to the left exits the system). When a toppling at L occurs, z(L) decreases by 1 and $z(L-1)$ increases by 1 (i.e. the particle that should have gone to the right stays at L).

For further reading on the Oslo model, we refer to [12, 37–41].

Assume that we start with a system containing many grains. Thanks to the Dirichlet boundary condition for h to the right, the system looses grains (h-particles) there until it reaches a critical slope. This critical slope is equivalent to a critical density of z-particles, which leave the system at the left ⁷ . We note here that sites with no z-particles or a single z-particle are always stable, sites with two z-particles could be stable or unstable, and particles with three or more z-particles are always unstable. This means that the critical slope must be between 1 and 2 (in fact, it equals approximately 1.8).

An important feature of the Oslo model is that it is Abelian (commutative). In our definition we explained which topplings are invoked, but we did not mention the order in which they occur. It is not difficult to see that the final configuration, after all topplings, does not depend on this order: z(i) is given by the number of incoming particles (the topplings at i − 1 and i + 1) and the number of outgoing particles (twice the topplings at i). Since a toppling at one site cannot render another site inactive, the total number of topplings at each site does not depend on the order, hence the final configuration also does not depend on the order.

Thanks to the Abelianity of the model, we are allowed to choose freely the order of topplings. To compare the Oslo model with the models presented above, we choose the same type of dynamics, making each site topple with rate 1 if $z(i)\gt z_{\mathrm{c}}(i)$.

The Oslo model has been studied extensively, and impressively accurate simulations of its critical exponents are known [12]. For the sake of this paper we wish to mention two of these exponents. First, the roughness exponent is conjectured [12] to be $\zeta = \frac{1}{4}$, that is, $\left \langle \left[h(i)-h(j)\right]^2 \right \rangle \sim |i-j|^{2\zeta}$. In the language of the z-particles, this is a hyperuniform state [42], i.e. the number of particles between i and j fluctuates as $|i-j|^\zeta \ll |i-j|^{1/2}$. The second relevant exponent is the dynamic exponent z_Oslo, conjectured [12] to be $z_\text{Oslo} = \frac{10}{7}$.

5.3.1. Local time and global time.

5.3.2. Hydrodynamic behaviour during the lifetime of an avalanche.

5.3.3. Comparing the Oslo model with aaEW.

As discussed in the previous section, we consider time scales much shorter than the lifetime of an avalanche. We can therefore assume that the system evolves according to a hydrodynamic limit of the type in equation (5), with an anchored boundary (at the left). In this model, the analogue of the particles in the DEP or the Toom model are the stress variables z.

Unlike the DEP or Toom's model, in the Oslo model not only the number of z-particles is conserved but also their center of mass, i.e.

$$\begin{equation} \partial_t \int_w^y h(x,t) \,\mathrm{d} x = \text{topplings at }w - \text{topplings at }y. \end{equation} \tag{ 54 }$$

Let us consider each of the terms in equation (5): The diffusion term $D \nabla^2 h$ and the advection term $\mu \nabla h$ integrate to

$$\begin{align} \int_w^y \mathrm{d} x \, D \nabla^2 h(x,t) &= D [\nabla h(y,t) - \nabla h(w,t)] , \end{align} \tag{ 55 }$$

$$\begin{align} \int_w^y \mathrm{d} x \, \mu\nabla h(x,t) &= \mu[ h(y,t) - h(w,t)], \end{align} \tag{ 56 }$$

both bounded uniformly in the interval length y − w (i.e. they remain bounded even if y − w becomes large). The quadratic and constant terms are given by

$$\begin{align} \int_w^y \mathrm{d} x \, \big[ \lambda_2 \nabla(h(x,t))^2 + c \big] \end{align} \tag{ 57 }$$

The zero-current condition (6) implies that this term vanishes on average. However its fluctuations grow with y − w, so it cannot be bounded uniformly. This means that $\lambda_2 = 0$. The cubic term

$$\begin{equation} \int_w^y \mathrm{d} x \, \lambda_3 [\nabla h(x,t)]^3 \end{equation} \tag{ 58 }$$

is also a fluctuating term which cannot be bounded uniformly, unless $\lambda_3 = 0$. In summary, the conservation of the center of mass forces $\partial_t h(x,t)$ to be a total derivatives, which is a stronger constraint than the particle-hole symmetry in Toom or the DEP ⁸ .

The argument above tells us that, for times shorter than $L^{z_\text{Oslo}}$, the Oslo model has an aaEW behaviour. Since $z_\text{Oslo} \gt z_\text{aaEW} = 1$, we are allowed to consider the model up to times $L^{z_\text{aaEW}}$. By then the surface reaches a roughness of $\zeta_\text{aaEW} = \frac{1}{4}$, which remains the roughness of the Oslo surface.

In figure 5 we show the auto-correlation function in the Oslo model, as a function of 'grain-time' t: the protocol is to add one grain, and then to let the full system relax during $n = 1/r$ iterations, or until all sites are stable. As a result, the time t in this simulation is the number of added grains, and r the driving rate. What we measure is the auto-correlation function in the steady state,

$$\begin{equation} \left\langle \left[ h(x,t)-h(x,0) \right]^2\right\rangle \equiv \left\langle \left[ h(x,t+t^{^{\prime}})-h(x,t^{^{\prime}}) \right]^2\right\rangle . \end{equation} \tag{ 59 }$$

One observes on figure 5 that for adiabatic driving (r → 0)

• (i)  
  For each x, it reaches a plateau after some time τ_x . The larger x, the larger τ_x .
• (ii)  
  A scaling collapse can be achieved by the ansatz

  $$\begin{equation} \left\langle \left[ h(x,t)-h(x,0) \right]^2\right\rangle \simeq x^{2\zeta} f_x\left(t/x^{z} \right) \quad \Longleftrightarrow \quad f_x(t) = \left\langle \left[ h(x,tx^z)-h(x,0) \right]^2\right\rangle, \end{equation} \tag{ 60 }$$

  where $f_x(t) \to f(t)$ for $x = 1,2,4,\ldots,L/2$.
• (iii)  
  The best scaling collapse is achieved by ζ = 0.258 (close to the predicted $\zeta = 1/4$), and z = 1.24, definitely smaller than $z = 10/7 = 1.428\,57$.
• (iv)  
  The last point x = L behaves differently, and its scaling function $f_L(t)$ does not collapse together with the others onto a master curve. While $f_x(t)\sim t^{2\zeta/z}$ with a power given by ζ = 0.258, and z = 1.24, the last point has a slope approaching $2\zeta/z = 7/20$.

Zoom In Zoom Out Reset image size

Let us now turn to a finite driving rate $r = 1/5$, i.e. after adding one grain on the left, we try to topple each site 5 times. Having a finite injection rate, t can both be interpreted as the number of injected grains, or time (divided by 5). We now observe that

• (i)  
  Using the theoretically predicted values of $\zeta = 1/4$ and $z = 10/7$ results in a decent scaling collapse, which improves for larger x.
• (ii)  
  The slope of f(t) in a logarithmic scale does not seem to approach $2\zeta/z$, except for the last point.
• (iii)  
  Things improve for larger system sizes.
• (iv)  
  Remarkably, while we studied smaller driving rates $r = 1/20$ and $r = 1/10$ (not shown) even at the large driving rate of $r = 1/5$ the plots are almost unchanged as compared to r = 0, indicating that we are still in a critical state.
• (v)  
  This critical state at large driving can be described by the anchored advected Edwards Wilkinson equation (17).

This confirms our findings that the roughness $\zeta = 1/4$ in the Oslo model has the same origin as in the anchored Edwards–Wilkinson equation (aaEW).

Here we give the relation to depinning in the Edwards–Wilkinson model, of a string of size L [2, 12, 22, 23]. The latter has an equation of motion,

$$\begin{equation} \partial_t u(i,t) = \nabla^2 u(i,t) + F(i,u(i,t)), \end{equation} \tag{ 61 }$$

which we read discretized in time t and space i ($\nabla ^2$ is the discrete Laplacian w.r.t. the variable i). The random forces $F(i,u)$ are uncorrelated Gaussian random variables with unit variance. To connect this to the Oslo sandpile, define

$$\begin{equation} u(i):= \sum_{j=i+1}^L h(i) + \# \{\textrm{grains fallen off at the right end}\} . \end{equation} \tag{ 62 }$$

As a consequence, the discrete Laplacian

$$\begin{equation} \nabla^2 u(i) = z(i), \end{equation} \tag{ 63 }$$

and the random force $F(i,u)$ stems from the fact that if $u(i)\to u(i)+1$, $F(i,u)\to F(i,u+1)$ is a new random variable, which identifies with $z_{\mathrm{c}}(i)$ in the Oslo model. The interpretation is that of a string pulled at site i = 0 (its left end), with Neumann boundary conditions (no force) at the right end. Its average profile is parabolic,

$$\begin{align} \left\langle u(i) \right\rangle \approx \frac{\left\langle z\right\rangle}{2} (L-i)^2 + \# \{\textrm{grains fallen off at the right} \}. \end{align} \tag{ 64 }$$

As the disorder is renewed after each displacement, it falls into the random-field universality class of depinning [2]. Since $\nabla u \sim h$, a roughness exponent $\zeta = \frac{1}{4}$ for h is the same as a roughness exponent

$$\begin{equation} \zeta^{\mathrm{qEW}}_{d=1} = \frac{5}{4} \end{equation} \tag{ 65 }$$

for u. This is what we wanted to show.