Anchored advected interfaces, Oslo model, and roughness at depinning

TL;DR

This paper investigates the roughness and dynamic exponents of various 1D advected systems with absorbing boundaries, revealing specific exponent values and connecting the Oslo model to elastic string depinning.

Contribution

It identifies the dynamic exponents for several advected models and links the Oslo sandpile to elastic string depinning, providing new insights into their universality classes.

Findings

The models exhibit roughness exponent ζ=1/4.

Dynamic exponents z=1, 2, 1/2 are observed depending on the model.

The Oslo model has a conjectured dynamic exponent z=10/7, implying ζ=5/4 for elastic string depinning.

Abstract

There is a plethora of 1-dimensional advected systems with an absorbing boundary: the Toom model of anchored interfaces, the directed exclusion process where in addition to diffusion particles and holes can jump over their right neighbor, simple diffusion with advection, and Oslo sandpiles. All these models share a roughness exponent of $\zeta=1/4$, while the dynamic exponent $z$ varies, depending on the observable. We show that for the first three models $z=1$, $z=2$, and $z=1/2$ are realized, depending on the observable. The Oslo model is apart with a conjectured dynamic exponent of $z=10/7$. Since the height in the latter is the gradient of the position of a disordered elastic string, this shows that $\zeta =5/4$ for a driven elastic string at depinning.

Equations

$$\begin{array}{rl}{\partial }_{t}h(x,t)& =D{\nabla }^{2}h(x,t)-\mu \nabla h(x,t)+\sigma \xi (x,t),\\ h(0,t)& =0,⟨\xi (x,t)\xi (x^{\prime },t^{\prime })⟩=\delta (x-x^{\prime })\delta (t-t^{\prime }).\end{array}$$

$$⟨h(x,t)h(x,0)⟩=\sqrt{\frac{x}{2\pi }}f(t/x^{1/2}).$$

$$\begin{array}{rl}{\partial }_{t}h(x,t)& =D{\nabla }^{2}h(x,t)-\mu \nabla h(x,t)+\sum _{n\ge 2}{\lambda }_n{\left[\nabla h(x,t)\right]}^2+\sigma \xi (x,t),\\ h(t,0)& =0,⟨\xi (x,t)\xi (x^{\prime },t^{\prime })⟩=\delta (x-x^{\prime })\delta (t-t^{\prime }).\end{array}$$

$$j(\rho )=2\left[{\rho }^2(1-\rho )-r(1-\rho {)}^2\rho \right]=2\rho (1-\rho )\left[\rho -r(1-\rho )\right].$$

$$j(\rho )=\rho (1-\rho )+2{\rho }^2(1-\rho )+3{\rho }^3(1-\rho )+\cdots -r\left[(1-\rho )\rho +2(1-\rho {)}^2\rho +\dots \right]=\frac{\rho }{1-\rho }-\frac{r(1-\rho )}{\rho }.$$

$$\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).$$

$$\begin{array}{rl}{\partial }_{t}h(x,t)=& D{\nabla }^{2}h(x,t)-\mu \nabla h(x,t)+\sigma \xi (x,t),\\ h(t,0)=& 0,⟨\xi (x,t)\xi (x^{\prime },t^{\prime })⟩=\delta (x-x^{\prime })\delta (t-t^{\prime }).\end{array}$$

$$\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<\xi(x,t)\xi(x^{\prime},t^{\prime})\right>=\delta(x-x^{\prime})\delta(t-t^{\prime}).\end{split}$$

$$c=-{\lambda }_2⟨[\nabla h{]}^2⟩.$$

$$⟨h(x,t+\mathrm{d}t)-h(x,t)⟩=(-p_1+p_2-p_3-p_4+p_5+p_6)\mathrm{d}t.$$

$$⟨h(x,t+\mathrm{d}t)-h(x,t)⟩$$

$$D$$

$$\mu$$

$${\lambda }_2$$

$${\sigma }^2\mathrm{d}t=⟨\frac{1}{k}{\left[\sum _{x=1}^{k}h(x,t+\mathrm{d}t)-h(x,t)\right]}^2⟩.$$

$$⟨\frac{1}{k}{\left(\sum _{x=1}^{k}h(x,t+\mathrm{d}t)-h(x,t)\right)}^2⟩=p_1+p_2+4p_3+4p_5.$$

$${\sigma }^2=2\rho (1-\rho )(1+2r+2\rho -2r\rho ).$$

$$D=\frac{4r+1}{r+1},\mu =\frac{2r}{r+1},{\lambda }_2=2(r-1),{\sigma }^2=\frac{2r(5r+1)}{(r+1{)}^3}.$$

$$D=\frac{5}{2},\mu =1,{\lambda }_2=0,{\sigma }^2=\frac{3}{2}.$$

$${\partial }_t\tilde{h}(x,t)$$

$$\begin{array}{rl}& {\partial }_{t}h(x,t)=D{\nabla }^{2}h(x,t)-\mu \nabla h(x,t)+\xi (x,t),h(x=0,t)=0.\\ & ⟨\xi (x,t)\xi (x^{\prime },t^{\prime })⟩=\delta (x-x^{\prime })\delta (t-t^{\prime }).\end{array}$$

$$h(x,t)={\int }_{y>0}{\int }_{t^{\prime }<t}{P}_{\mu }(x,t|y,t^{\prime })\xi (y,t^{\prime }).$$

$$P_0(x,t)$$

$$P_0^{\mu }(x,t)$$

$$P^{\mu }(x,t|y,t^{\prime })$$

$$⟨h(x,t{)}^2⟩$$

$$P^{\mu }(x,t|y,0)$$

$$\left<h(x)^{2}\right>_{\mu}=\frac{\sigma^{2}}{\mu}\left<h\Big{(}\frac{\mu x}{D},\frac{\mu^{2}t}{D}\Big{)}^{2}\right>_{\mu=D=\sigma=1}.$$

$${\int }_{y>0}{\int }_{t>0}{P}^{\mu =1}(x,t|y,0{)}^2={\int }_{k,p}\frac{-16kp{\mathrm{e}}^{x[1+i(k+p)]}}{[2(k^2+p^2)+1][(k-p{)}^2+1][(k+p{)}^2+1]}.$$

$${⟨h(x{)}^2⟩}_{\mu =1}={\int }_{-\infty }^{\infty }\frac{\mathrm{d}p}{2\pi i}\frac{16p{\mathrm{e}}^{-\sqrt{p^2+\frac{1}{2}}x+ipx+x}}{(4p^2+1{)}^2}.$$

$${\int }_0^{\infty }\mathrm{d}tt{\mathrm{e}}^{-at}=\frac{1}{a^2},$$

$${\int }_{s>0}\frac{{\mathrm{e}}^{-as-\frac{1}{4s}}}{2\sqrt{\pi }{s}^{3/2}}={\mathrm{e}}^{-\sqrt{a}}.$$

$${⟨h(x{)}^2⟩}_{\mu =1}={\int }_{s>0}{\int }_{t>0}\frac{2tx{\mathrm{e}}^{-\frac{x^2}{4s{x}^2+16t}-\frac{s{x}^2}{2}-\frac{1}{4s}-t+x}}{\pi s^{3/2}(s{x}^2+4t{)}^{3/2}}.$$

Full text

1 Introduction

Interfaces subject to quenched disorder describe a variety of physical phenomena [1, 2], such as domain walls in magnets [3, 4, 5], contact-line depinning [6], fracture [7, 8], or earthquakes [9]. Two universality classes have to be distinguished: equilibrium and depinning. Though many numerical studies exist [10, 11, 12], and a field theory was developed [13, 14, 2], there are few exact results.

A notable exception in equilibrium is the roughness exponent $\zeta_{\rm RB}^{d=1}=2/3$ for a $1+1$ dimensional directed polymer in random-bond (RB) disorder, itself related to the KPZ universality class [15, 1] with roughness $1/2$ and dynamic exponent $z=1/\zeta_{\rm RB}^{d=1}=3/2$ [16]. For random-field (RF) disorder of a $d$-dimensional interface in $d+1$ dimensions in equilibrium, scaling arguments correctly predict $\zeta_{\rm RF}^{d}=(4-d)/3$, which can experimentally be seen even in dimension $d=0$ [17], where it reduces to the Sinai model [18].

At depinning, there is the single random-field universality class, and no analytic result is known, apart from $d=0$ for which $\zeta_{\rm dep}^{d=0}=2^{-}$ [19]. Based on numerical simulations, it was recently conjectured [12] that a driven 1-dimensional string has a roughness exponent of $\zeta_{\rm dep}^{d=1}=5/4$. Here we aim at demonstrating this value. A key observation is that if the roughness of the string at depinning is $\zeta_{\rm dep}^{d=1}=5/4$, the Oslo model [20, 21], of which the height can be viewed as the gradient of the position of a driven string pulled at one end, has a roughness of $\zeta_{\rm Oslo}^{d=1}=1/4$.

There are several 1-dimensional systems with a scaling exponent of $1/4$, but this is generally in the temporal domain. E.g. a marked monomer on a polymer grows with time as $\left<[h(0,t)-h(0,t^{\prime})]^{2}\right>\sim|t-t^{\prime}|^{2H}$, with $H=1/4$. This is obtained from standard scaling arguments as $H=\zeta/z$, with $\zeta=1/2$ and $z=2$. If these systems are advected, i.e. $z=1$, the roughness exponent becomes $\zeta=Hz\to 1/4$ [22, 23, 24, 25, 26].

Here we consider a class of models which allow for a hydrodynamic description equivalent to diffusion of a scalar field $h(x,t)$ combined with advection away from an absorbing boundary at $x=0$:

[TABLE]

We call it the anchored advected Edwards-Wilkinson equation (aaEW). While in the bulk the stationary state $h(x,\infty)$ can be obtained from a Brownian motion in the comoving frame, close to the boundary we establish in section 3.5 that $\left<h(x,\infty)^{2}\right>\sim\sqrt{x}$, thus $\zeta=1/4$.

This system possesses three distinct dynamical scaling behaviours, characterized by a dynamical exponent $z=1$ due to the advection, $z=2$ for bulk properties in the comoving frame, and $z=1/2$ for the decay of the equal-point correlation function (in the steady state),

[TABLE]

The function $f$ is obtained analytically in section 3.7.

The model (aaEW) is the proper thermodynamic description for a variety of systems. We give two examples – the directed exclusion process (section 2.1) and the Toom interface model (section 2.2). For the Oslo model (section 5) a similar mapping can be constructed, with a crucial difference: The local time which sets the clock for an update is advanced by the toppling itself, leading to a different dynamical exponent, conjectured to be $z=10/7$ [12], thus larger than the corresponding exponent $z=1$ in aaEW, but smaller than the one for diffusion ($z=2$). Still, the roughness exponent equals $\zeta=1/4$, independent of whether one considers bulk or boundary driving. The Oslo model (section 5) achieves this by propagating out from the seed in both directions; the seed position effectively acts as an absorbing boundary.

In the derivation of a hydrodynamic description as in Eq. (aaEW), additional terms may appear,

[TABLE]

The perturbation with $\lambda_{2}$ is the most relevant one, leading to the anchored advected Kardar-Parisi-Zhang universality class (aaKPZ) discussed in section 4. If this term is forbidden by symmetry, the next relevant one is $\lambda_{3}$. It is marginally relevant, and believed to modify the effective parameters in the symmetric case [25], but it cannot change the stationary measure in the bulk (section 4.2).

This paper is organized as follows: In section 2 we introduce models in the aaEW and aaKPZ universality classes, first the directed exclusion process (section 2.1), then Toom’s interface model (section 2.2); their hydrodynamic description is obtained in section 2.3. In section 3 we derive analytic results for the aaEW universality class, followed by results for the aaKPZ class in section 4. The Oslo model is treated in section 5.

2 Models in the aaEW and aaKPZ universality classes

2.1 The Directed Exclusion Process (DEP)

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

10 $\to$ 01 with rate $1$, 2. 2.

01 $\to$ 10 with rate $1$, 3. 3.

110 $\to$ 011 with rate $1$, 4. 4.

001 $\to$ 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\neq 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

[TABLE]

Therefore, $\rho_{\rm 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_{\rm c}=\frac{r}{1+r}$ in the domain of observation.

2.2 Toom’s interface model

Toom’s model is the subject of many theoretical and numerical studies [27, 22, 28, 29, 30]. 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 flips 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\neq 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, and the total particle current is

[TABLE]

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

2.3 Continuum theory

As suggested in [22, 29], 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

[TABLE]

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

[TABLE]

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

[TABLE]

Since $\left<h(x,t)\right>=0$ by construction, the first two terms on the r.h.s. vanish on average, and hence

[TABLE]

An additional term with $\lambda_{3}\neq 0$ may lead to logarithmic corrections [31].

A direct derivation of this limit for Toom’s model appears in [22, 29]. We describe here the derivation for the DEP.

2.4 Continuum limit of 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 [math] otherwise. Then the height function $h$ is defined as $h(x,t):=\sum_{y=1}^{x}[\eta(y,t)-\rho]$.

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

$h(x,t+{\rm d}t)=h(x,t)-1$ with probability $p_{1}=\eta(x,t)[1-\eta(t,x+1)]{\rm d}t$, 2. 2.

$h(x,t+{\rm d}t)=h(x,t)+1$ w.p. $p_{2}=[1-\eta(x,t)]\eta(t,x+1){\rm d}t$, 3. 3.

$h(x,t+{\rm d}t)=h(x,t)-1$ w.p. $p_{3}=\eta(x,t)\eta(x+1,t)[1-\eta(x+2,t)]{\rm d}t$, 4. 4.

$h(x,t+{\rm d}t)=h(x,t)-1$ w.p. $p_{4}=\eta(x-1,t)\eta(x,t)[1-\eta(x+1,t)]{\rm d}t$, 5. 5.

$h(x,t+{\rm d}t)=h(x,t)+1$ w.p. $p_{5}=r[1-\eta(x,t)][1-\eta(x+1,t)]\eta(x+2,t){\rm d}t$, 6. 6.

$h(x,t+{\rm d}t)=h(x,t)+1$ w.p. $p_{6}=r[1-\eta(x-1,t)][1-\eta(x,t)]\eta(x+1,t){\rm d}t$.

First, we calculate

[TABLE]

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

[TABLE]

The noise term is obtained from

[TABLE]

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

[TABLE]

Note that $p_{4},p_{6}$ are spatial shifts of $p_{3},p_{5}$ hence we must count them only once. This yields

[TABLE]

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

[TABLE]

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

[TABLE]

3 Particle-hole symmetric case – the aaEW universality class

3.1 Periodic BC versus anchored interface

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. 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,

[TABLE]

(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.

3.2 Bulk behaviour

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 [32, 30], and we expect a similar behaviour in all models of this universality class.

3.3 Edge behaviour

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 [22, 29, 26]. The combination of the Dirichlet boundary and the advection term transfers temporal correlations into spacial ones. In the following, we derive analytic results for the edge behaviour in the anchored case.

3.4 Basic formulas

We solve the advected diffusion equation with Dirichlet boundary conditions, see also [26]111In [26] $\mu$ has opposite sign, which is equivalent to studying the equation on the negative half line $(-\infty,0]$.. By simple rescaling we write equation (aaEW) in the form:

[TABLE]

Eq. (17) for $x>0$ is solved by

[TABLE]

Here

[TABLE]

3.5 Analytic result for the variance of the height

The roughness exponent $\zeta$ is given by the behaviour of $\left<h(x,t)^{2}\right>$ in the steady state,

[TABLE]

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

[TABLE]

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

[TABLE]

Then

[TABLE]

Integrating Eq. (23) over $k$ yields

[TABLE]

Eq. (24) can be simplified, using for the denominator

[TABLE]

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

[TABLE]

Then the integral over $p$ can be done, leading to

[TABLE]

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

[TABLE]

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

[TABLE]

The integral is known analytically,

[TABLE]

The function $\boldsymbol{L}_{n}(x)$ is the modified Struve-function. For $\mu\neq 1$ this reads according to Eq. (22)

[TABLE]

One can approximate the Bessel function for large argument as

[TABLE]

Then the integration can be done analytically,

[TABLE]

This gives a roughness exponent

[TABLE]

3.6 Approximate calculation for large $x$

From Figure 1 we see that, for large $x$, we can approximate Eq. (18) as

[TABLE]

Then

[TABLE]

This reproduces the result of Eq. (33).

3.7 Decay of the auto-correlation function

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

[TABLE]

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

[TABLE]

Therefore

[TABLE]

The $\Gamma$-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,

[TABLE]

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$

[TABLE]

Putting back the dependence on $\mu$, $D$ and $\sigma$ yields with Eq. (22)

[TABLE]

3.8 Dynamic exponents $z=2$, $z=1$ and $z=1/2$

Our calculations above teach us that

[TABLE]

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

3.9 Numerical verification

In the symmetric case, we expect the (centered) 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 Eq. (42).

In Fig. 2 we show this comparison for Toom’s model. It is worth noting, as already mentioned in [22], 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 $\mu$, $D$, and $\sigma$). For more thorough numerical studies of this model we refer the reader to [22, 29, 33].)

In Fig. 3 we show the same plot for the symmetric DEP. There are noticeable finite-size effects, but for large systems the rescaled correlator $\left<h(x,t)h(0,0\right>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 Eq. (17) to account for the coefficients $D$ and $\sigma$ in Eq. (aaEW), and plugging in $D,\mu$ and $\sigma$ found in Section 2.4, we obtain

[TABLE]

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

[TABLE]

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 Refs. [31, 33]. 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.

The code used to simulate the DEP is available as a part of the arXiv version of this paper.

3.10 Discrete argument

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 $\alpha>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 $\alpha$ scale as $L^{-1/2}$, and hence $z=\frac{1}{2}$.

4 Asymmetric case – the aaKPZ universality class

4.1 Introduction

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 [22, 23, 24, 25], 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$.

4.2 Protection of the stationary measure

Consider the equation of motion (17), with a subleading KPZ-term, and a subsubleading term $\sim(\nabla h)^{3}$, 222Here we extend an argument initially given for the KPZ equation. While it may be well known in this more general form, we could not find a reference in the literature. J. Krug mentions it in his Diplomarbeit of 1985 at Munich university.

[TABLE]

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

[TABLE]

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}^{\rm ss}[h]=0$ can be found by asking that

[TABLE]

This is solved by

[TABLE]

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

[TABLE]

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 $\mu$ and $\lambda_{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 $\eta$, $D$, $\mu$, and $\lambda_{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<h(x,t)h(x,0)\right>$ depends on the distance $x$ to the boundary, they could appear there. Since the static 2-point function $\left<h(x,t)^{2}\right>$ 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 rigorous.

5 Oslo model

5.1 Definition

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

[TABLE]

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

[TABLE]

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$,

[TABLE]

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

[TABLE]

In its original version, the random number is $1$ or $2$ with probability $1/2$. To obtain Fig. 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)>z_{\rm 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)>z_{\rm 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 Refs. [34, 35, 12, 36, 37, 38].

5.2 Phenomenology of the Oslo model

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.333Note that since due to Eq. (50) $z$ is the derivative of $h$. Thus Neumann boundary conditions for $h$ are Dirichlet boundary conditions for $z$, and vice versa. 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. In order 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)>z_{\rm 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 [39], 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_{\text{Oslo}}$, conjectured [12] to be $z_{\text{Oslo}}=\frac{10}{7}$.

5.3 Why is $\zeta=\frac{1}{4}$ in the Oslo model?

5.3.1 Local time and global time

The quenched noise of the Oslo model may be thought of as a local time shift: $z_{c}$ is only updated after a toppling, while an annealed version would update $z_{c}$ at a fixed rate. We therefore define the local time at a site $x$ as the number of toppling events at $x$. This delays the topplings as compared to the models discussed in sections 3, but when the topplings occurr, their dynamics is faster, so as to “catch up”. As a result, we expect the dynamical critial exponent $z$ to be smaller than $2$, and indeed $z_{\text{Oslo}}\approx\frac{10}{7}<2$.

5.3.2 Hydrodynamic behaviour during the lifetime of an avalanche

A key observation for the Oslo model is that, during an avalanche, the local time and the global time propagate approximately at the same rate. Therefore, as long as we are interested in the evolution of an avalanche at times which are much shorter than its survival time, we may replace the quenched noise with an annealed one.

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 directed exclusion process444We already point out that this is related to the fact that $h(x,t)$ can be thought of as the gradient of a field $u(x,t)$, see Eq. (62)..

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}}>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.

5.4 Adiabatic versus finite advection in the Oslo model

In Fig. 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,

[TABLE]

One observes on Fig. 5 that for adiabatic driving ($r\to 0$)

1. (i)

For each $x$, it reaches a plateau after some time $\tau_{x}$. The larger $x$, the larger $\tau_{x}$. 2. (ii)

A scaling collapse can be achieved by the ansatz

[TABLE]

where $f_{x}(t)\to f(t)$ for $x=1,2,4,...,L/2$. 3. (iii)

The best scaling collapse is achieved by $\zeta=0.258$ (close to the predicted $\zeta=1/4$), and $z=1.24$, definitely smaller than $z=10/7=1.42857$. 4. (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 $\zeta=0.258$, and $z=1.24$, the last point has a slope approaching $2\zeta/z=7/20$.

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

1. (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$. 2. (ii)

The slope of $f(t)$ in a logarithmic scale does not seem to approach $2\zeta/z$, except for the last point. 3. (iii)

Things improve for larger system sizes. 4. (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. 5. (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).

5.5 Oslo and quenched Edwards-Wilkinson

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

[TABLE]

which we read discretized in time $t$ and space $u$. The random forces $F(i,u)$ are uncorrelated Gaussian random variables with unit variance. To connect this to the Oslo sandpile, define

[TABLE]

As a consequence, the discrete Laplacian

[TABLE]

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_{\rm 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,

[TABLE]

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

[TABLE]

for $u$. This is what we wanted to show.

6 Conclusion

As we have shown, there is a large class of models which have a roughness exponent of $\zeta=1/4$. This encompasses the Toom model, advected diffusion with an absorbing boundary, the symmetric directed exclusion process, and the Oslo model. Observing that the height in the Oslo model corresponds to the slope of the height in the quenched Edwards-Wilkinson model, we showed that the roughness exponent of the latter is $\zeta_{\rm qEW}=5/4$. This puts onto a firm basis the recent conjecture [12] that the numerically observed exponent of $\zeta_{d=1}^{\rm qEW}=1.25$ [11] is exactly $\zeta_{d=1}^{\rm qEW}=5/4$.

While these models have the same roughness exponent, their dynamical exponent $z$ can be quite different, and depend on the observable. We identified $z$ for the first three models, the anchored advected universality class, where we observed a dynamical exponent of $z=2$, $z=1$ and even $z=1/2$, depending on the obversable. Our correspondence to the Oslo model cannot (yet) provide its dynamic exponent $z$, conjectured to be $z=10/7$, since the arguments use both a global and a local time, which are distinct in this case. Given that $z=10/7$ is again a simple fraction, we hope that an exact argument may be found as well.

Acknowledgements

We are grateful to Joachim Krug for explaining how advection can convert temporal correlations into spatial ones, and for directing us to the relevant literature. We thank Peter Grassberger for a useful exchange.