L-6

Goal: We solve the mean field version of the cellular automaton, derive its avalanche statistics and make a connection with the Bienaymé-Galton-Watson process used to describe an epidemic outbreak.

Fully connected (mean field) model for the cellular automaton

Let's study the mean field version of the cellular automata introduced in the previous lecture. We introduce two approximations:

• Replace the Laplacian, which is short range, with a mean field fully connected interction

$\sigma _{i}=h_{CM}-h_{i}+m^{2}(w-h_{i}),\quad$

.

• The local threshold are all equal. In particular we set

$\sigma _{i}^{th}=1,\quad \forall i$

.

As a consequence, in the limit $L\to \infty$, the statistical properties of the system are described by the distribution of the local stresses $\sigma _{i}$. For simplicity, instead of the stresses, we study the distance from threshold

$x_{i}=1-\sigma _{i}$

The instability occurs when a block is at $x_{i}=0$ and is followed by its stabilization and a redistribution on all the blocks :

${\begin{cases}x_{i}=0\to x_{i}=\Delta &(stabilization)\\x_{j}\to x_{j}-{\frac {1}{L}}{\frac {\Delta }{1+m^{2}}}&(redistribution)\\\end{cases}}$

Dynamics

Our goal is thus to determine the distribution $P_{w}(x)$ of all blocks, given their intial distribution, $P_{0}(x)$, and a value of $w$. Let's decompose in steps the dynamics

• Drive: Increasing $w\to w+dw$ each block decreases its distance to threshold

$x_{i}\to x_{i}-m^{2}dw$

.

As a consequence

$P_{w+dw}(x)=P_{w}(x+m^{2}dw)\sim P_{w}(x)+m^{2}dw\partial _{x}P_{w}(x)$

• Stabilization : A fraction $m^{2}dwP_{w}(0)$ of the blocks is unstable. The stabilization induces the change $m^{2}dwP_{w}(0)\to m^{2}dwP_{w}(0)g(x)$. Hence, one writes

$\partial _{w}P_{w}(x)\sim m^{2}\left[\partial _{x}P_{w}(x)+P_{w}(0)g(x)\right]$

The stabilization of the unstable blocks induce a drop of the force per unit length

$m^{2}dwP_{w}(0)\int dxxg(x)=m^{2}dwP_{w}(0){\overline {\Delta }}$

\

• Redistribution This drop is (partially) compensated by the redistribution. The force acting on all points is increased:

$x_{i}\to x_{i}-m^{2}dwP_{w}(0){\frac {\overline {\Delta }}{1+m^{2}}}$

Again, most of the distribution will be driven to instability while a fraction of the blocks become unstable... we can write

$\partial _{w}P_{w}(x)=m^{2}\left[\partial _{x}P_{w}(x)+P_{w}(0)g(x)\right]\left[1+P_{w}(0){\frac {\overline {\Delta }}{1+m^{2}}}+(P_{w}(0){\frac {\overline {\Delta }}{1+m^{2}}})^{2}+\ldots \right]$

and finally:

$\partial _{w}P_{w}(x)={\frac {m^{2}}{1-P_{w}(0){\frac {\overline {\Delta }}{1+m^{2}}}}}\left[\partial _{x}P_{w}(x)+P_{w}(0)g(x)\right]$

Stationary solution

Increasing the drive the distribution converge to the fixed point:

$0=\partial _{x}P_{\text{stat}}(x)+P_{\text{stat}}(0)g(x)$

• Determine $P_{\text{stat}}(0)={\frac {1}{\overline {\Delta }}}$ using

$1=\int _{0}^{\infty }dx\,P_{\text{stat}}(x)=-\int _{0}^{\infty }dx\,x\partial _{x}P_{\text{stat}}(x)$

• Show

$P_{\text{stat}}(x)={\frac {1}{\overline {\Delta }}}\int _{x}^{\infty }g(z)dz$

which is well normalized.

Critical Force

The average distance from the threshold gives a simple relation for the critical force, namely $1-f_{c}={\overline {x}}$. Hence for the automata model we obtain:

$f_{c}=1-\int _{0}^{\infty }dxxP_{\text{stat}}(x)=1-{\frac {1}{2}}{\frac {\overline {\Delta ^{2}}}{\overline {\Delta }}}$

Exercise:

Let's assume an exponential distribution of the thresholds and show

• $P_{\text{stat}}(x)=e^{-x/{\overline {\Delta }}}/{\overline {\Delta }}$
• $f_{c}=1-{\overline {\Delta }}$

Avalanches or instability?

We consider an avalanche starting from a single unstable site $x_{0}=0$ and the sequence of sites more close to instabitity $x_{1}<x_{2}<x_{3}\ldots$. For each unstable block, all the blocks receive a random kick:

${\frac {\Delta _{1}}{(1+m^{2})L}},\quad {\frac {\Delta _{2}}{(1+m^{2})L}},\quad {\frac {\Delta _{3}}{(1+m^{2})L}},\ldots$

with $\Delta _{1},\Delta _{2},\Delta _{3},\ldots$ drwan from $g(\Delta )$ Are these kick able to destabilize other blocks?

Given the initial condition and $w$, the state of the system is described by $P_{w}(x)$. From the extreme values theory we know the equation setting the average position of the most unstable block is

$\int _{0}^{x_{1}}P_{w}(t)dt={\frac {1}{L}}$

Hence, for large systems we have

$x_{1}\sim {\frac {1}{LP_{w}(0)}},\;x_{2}\sim {\frac {2}{LP_{w}(0)}},\;x_{3}\sim {\frac {3}{LP_{w}(0)}},\ldots$

Hence we need to compare the mean value of the kick with the mean gap between nearest unstable sites:

${\frac {\overline {\Delta }}{(1+m^{2})L}}\quad {\text{versus }}\quad {\frac {1}{P_{w}(0)L}}$

Note that $L$ simplifies. We expect three possibilities:

• if the mean kick is smaller than the mean gap the system is subcritical and avalanches quickly stops.
• if the mean kick is equal to the mean gap the system is critical and avalanches are power law distributed
• if the mean kick is larger of the mean gap the system is super-critical and avalanches are unstable.

Note that in the stationary regime the ratio between mean kick and mean gap is $1/(1+m^{2})$. Hence, the system is subcritical when $m>0$ and critical for $m=0$

Mapping to the Brownian motion

Let's define the random jumps and the associated random walk

$\eta _{1}={\frac {\Delta _{1}}{(1+m^{2})L}}-x_{1},\;\eta _{2}={\frac {\Delta _{2}}{(1+m^{2})L}}-(x_{2}-x_{1}),\;\eta _{3}={\frac {\Delta _{3}}{(1+m^{2})L}}-(x_{3}-x_{2})\ldots$ $X_{n}=\sum _{i=1}^{n}\eta _{i}\quad \quad {\text{with}}\;{\overline {\eta _{i}}}={\frac {\overline {\Delta }}{L(1+m^{2})}}-{\frac {1}{LP_{w}(0)}}$

An avalanche is active until $X_{n}$ is positive. Hence, the size of the avalanche identifies with first passage time of the random walk.

• Critical case : In this case the jump distribution is symmetric and we can set $X_{0}=0$. Under these hypothesis the Sparre-Andersen theorem state that the probability that the random walk remains positive for $n$ steps is independent on the jump disribution and for a large number of steps becomes $Q(n)\sim {\frac {1}{\sqrt {\pi n}}}$. Hence, the distribution avalanche size is

$P(S)=Q(S)-Q(S+1)\sim {\frac {1}{\sqrt {\pi S}}}-{\frac {1}{\sqrt {\pi (S+1)}}}\sim {\frac {1}{2{\sqrt {\pi }}}}{\frac {1}{S^{3/2}}}$

This power law is of Gutenberg–Richter type. The universal exponent is $\tau =3/2$

• Stationary regime: Replacing ${\frac {1}{LP_{w}(0)}}$ with ${\frac {1}{LP_{\text{stat}}(0)}}={\frac {\overline {\Delta }}{L}}$ we get $\;{\overline {\eta _{i}}}\sim -{\frac {m^{2}}{1+m^{2}}}{\frac {\overline {\Delta }}{L}}$. For small m, the random walk is only sliglty tilted. The avalanche distribution will be power law distributed with $\tau =3/2$ until a cut-off

$S_{\max }\sim m^{-4}$

Bienaymé Galton Watson process

A time $t=0$ appears as infected individual which dies with a rate $a$ and branches with a rate $b$. On average, each infection generates in average $R_{0}=b/a$ new ones. Real epidemics corresponds to $R_{0}>1$.

At time $t$, the infected population is $n(t)$, while the total infected population is

$N(t)=\int _{0}^{t}n(t')dt'$

Our goal is to compute $P(N(t))$ and we introduce its Laplace Transform:

$Q_{s}(t)=\int _{0}^{\infty }P(N)e^{-sN}dN=\left\langle e^{-s\int _{0}^{t}n(t')dt'}\right\rangle$

. Note that the normalization imposes $Q_{0}(t)=1$.

• Evolution equation: Consider the evolution up to the time $t+dt$ as a first evolution from $0$ to $dt$ and a following evolution from $dt$ to $t+dt$. Derive the following equation for $Q_{s}(t)$

$Q_{s}(t+dt)=(1-(a+b)dt)e^{-sdt}Q_{s}(t)+adt+bdtQ_{s}^{2}(t)+O(dt^{2})$

which gives

${\frac {dQ_{s}(t)}{dt}}=-(a+b+s)Q_{s}(t)+a+bQ_{s}^{2}(t)$

• Critical case: the stationary solution: Let's set $b=a$ and $a=1$ to recover the results of the mean field cellular automata. In the limit $t\to \infty$ the total population coincides with the avalanche size, $N(t\to \infty )=S$. The Laplace transform of $P(S)$ is

$0=-(2+s)Q_{s}^{\text{stat}}+1+(Q_{s}^{\text{stat}})^{2}$

which gives

$Q_{s}^{\text{stat}}={\frac {(2+s)-{\sqrt {s^{2}+4s}}}{2}}\sim 1-{\sqrt {s}}+O(s)$

with

$\int _{0}^{\infty }dSP(S)e^{-sS}=Q_{s}^{\text{stat}}$

• Critical case: Asymptotics: We want to predict the power law tail of the avalanche distribution $P(S)\sim A\cdot S^{-\tau }$. Taking the derivative with respect to $s$ we have

$A\int _{0}^{\infty }dSS^{1-\tau }e^{-sS}={\frac {1}{2{\sqrt {s}}}}$

and conclude that $\tau =3/2$ and

$A={\frac {1}{2\int _{0}^{\infty }dze^{-z}/{\sqrt {z}}}}={\frac {1}{2{\sqrt {\pi }}}}$

Hence we find back our previous result

$P(S)\sim {\frac {1}{2{\sqrt {\pi }}}}{\frac {1}{S^{3/2}}}$