2. Basic estimates for a subcritical network

The aim of this section is to give bounds on the expected number of excited nodes at a given time, the total number of nodes excited during the entire duration of the avalanche (referred to in what follows as the cumulative size of the avalanche), and the probability that a given node is excited at a fixed time $k\in\mathbb{N}.$ Throughout the section we consider a single network; that is, we assume that $n\in\mathbb{N}$ and $p\in (0,1)$ are fixed in (1.2). Most of the results here are of an auxiliary nature, but some, in particular Propositions 2.2, 2.3, and 2.4, appear to be of independent interest. Proposition 2.2 is concerned with the evolution of a measure of heterogeneity for the network; Proposition 2.3 gives an upper bound, uniform on $V_n$, on the probability that a given node $x\in V_n$ is excited at time k; and Proposition 2 gives tight lower and upper bounds for the expected cumulative size of the avalanche in the subcritical regime. For future convenience, we formulate our results in terms of arbitrary bounds $c>0$ and $d>0$ that satisfy the following condition.

Condition 2.1 $0<d\leq pn \leq c.$

The next series of propositions is formulated for an arbitrary $c>0,$ even though in the sequel the results are used primarily in the case when $c\in (0,1).$ As we will see in the next section, the case $c\in (0,1)$ corresponds to a subcritical regime of the avalanche model, which turns out to be p satisfying the condition $pn<1.$

Let $(\mathcal{F}_k)_{k\in\mathbb{Z}_+}$ be the natural filtration for the sequence $(X_k)_{k\in\mathbb{Z}_+}$; that is, $\mathcal{F}_k$ is the $\sigma$-algebra generated by $X_0,\ldots,X_k.$ It follows from (1.2) that with probability one,

(2.1) \begin{eqnarray}\mathbb{E}(X_{k+1}\,|\,\mathcal{F}_k)\leq (n-X_k)\Bigl[1-\Bigl(1-\frac{c}{n}\Bigr)^{X_k}\Bigr]\leq cX_k\frac{n-X_k}{n}\leq cX_k.\end{eqnarray}

Notice that the formula remains formally true when $X_k=0.$ Thus, we have the following.

We remark that for $c\in [1,4),$ the first inequality in (2.1) produces an alternative sensible estimate $\mathbb{E}(X_k)\leq \frac{cn}{4}$ for all $k\in\mathbb{N}.$ These straightforward estimates will be improved for the entire range of the network parameters (i.e. for np less than, greater than, or equal to 1) in Section 6 below.

For $k\in\mathbb{Z}_+,$ let

(2.2) \begin{eqnarray}H_k=X_k(n-X_k).\end{eqnarray}

The following is immediate from Proposition 2.1 and (2.1).

Corollary 2.1. Let Condition 2.1 hold. Then

\begin{equation*}\mathbb{E}(X_m)\leq \frac{c^m}{n}\mathbb{E}(H_0)\leq c^m\mathbb{E}(X_0)\end{equation*}

for all $m\in\mathbb{N}.$

Indeed, Proposition 2.1 implies that $\mathbb{E}(X_m)\leq c^{m-1}\mathbb{E}(X_1),$ and a subsequent application of (2.1) with $k=0$ yields the corollary.

For $k\in\mathbb{Z}_+$ and $x\in V_n,$ let

(2.3) \begin{eqnarray}\mathcal{E}_k(x)=\left\{\begin{array}{l@{\quad}l}1&\textrm{if}\ x\in A_k,\\ \\[-7pt] 0&\textrm{if}\ x\not\in A_k,\end{array}\right.\end{eqnarray}

be the indicator of the event that a given node $x\in V_n$ is excited at time $k\in\mathbb{Z}_+.$ Note that $H_k=X_k(n-X_k)$ can be interpreted as a (non-normalized) measure of the heterogeneity of the network, because

\begin{eqnarray*}H_k&=&X_k(n-X_k)=\Bigl(\sum_{x \in V_n} \mathcal{E}_k(x)\Bigr)\cdot \Bigl(\sum_{x \in V_n}\bigl(1-\mathcal{E}_k(x)\bigr)\Bigr)\\&=&\sum_{x,y \in V_n} \mathcal{E}_k(x)\bigl(1-\mathcal{E}_k(y)\bigr),\end{eqnarray*}

and thus

\begin{equation*}\frac{2H_k}{n(n-1)}\end{equation*}

is the probability that two nodes randomly chosen at time k are not in the same state. As discussed in Section 6, it turns out that in the subcritical case $np\in(0,1)$ as well as in the critical case $np=1,$ for large values of n, the number of excited nodes loosely speaking decreases to zero almost monotonically. This observation motivates the following result, which in particular implies that the expected heterogeneity decreases to zero monotonically in the subcritical regime (see also Corollary 6.2 below).

Proof. It follows from (1.2) that

(2.4) \begin{eqnarray}\mathbb{E}\big(X_{k+1}^2|X_k\big)&=&\big(n-X_k\big)\big(1-q^{X_k}\big)q^{X_k}+\big(n-X_k\big)^2\big(1-q^{X_k}\big)^2.\end{eqnarray}

Notice that the formula remains true when $X_k=0.$

It follows that

(2.5) \begin{eqnarray}\nonumber&&\mathbb{E}\bigl[X_{k+1}(n-X_{k+1})\bigr|X_k\bigr]=n\mathbb{E}(X_{k+1}|X_k)-\mathbb{E}\big(X_{k+1}^2|X_k\big)\\[2pt]\nonumber &&\qquad=n(n-X_k)\big(1-q^{X_k}\big)- (n-X_k)\big(1-q^{X_k}\big)q^{X_k}-(n-X_k)^2\big(1-q^{X_k}\big)^2\\[2pt]\nonumber &&\qquad=(n-X_k)\big(1-q^{X_k}\big)\bigl[n-q^{X_k}-(n-X_k)\big(1-q^{X_k}\big)\bigr]\\[2pt]\nonumber &&\qquad=(n-X_k)\big(1-q^{X_k}\big)\bigl[X_k(1-q^{X_k})+(n-1)q^{X_k}\bigr]\\[2pt]\nonumber &&\qquad\leq(n-1)(n-X_k)\big(1-q^{X_k}\big)\leq n(n-X_k)\bigg(1-\left(1-\frac{c}{n}\right)^{X_k}\bigg)\\[2pt] &&\qquad\leq cX_k(n-X_k),\end{eqnarray}

where in the first inequality we used the fact that $X_k\leq n-1,$ and hence

\begin{eqnarray*}X_k\big(1-q^{X_k}\big)+(n-1)q^{X_k}\leq n-1.\end{eqnarray*}

The proof of the proposition is complete.

Recall $\mathcal{E}_k(x)$ from (2.3). Let

(2.6) \begin{eqnarray}\xi_k(x)= \mathbb{P}\bigl(\mathcal{E}_k(x)=1\bigr), \qquad x\in V_n,\end{eqnarray}

and

\begin{eqnarray*}\eta_k=\mathbb{P}(X_k>0)=\mathbb{P}(T>k), \qquad k \in\mathbb{Z}_+,\end{eqnarray*}

where T is the duration of the avalanche defined in (1.3). We have the following:

Proposition 2.3. Let Condition 2.1 hold. Then

\begin{equation*}\xi_k(x)\leq \frac{c^k}{n^2}\mathbb{E}(H_0)\end{equation*}

for all $k\in\mathbb{N}$ and $x\in V_n.$

The proof of the proposition is given below in this section, after the statement of a corollary. The claim is trivial if the distribution of $X_0$ is invariant with respect to permutation of nodes (see the proof below), but takes slightly more effort to establish when the symmetry is broken and the nodes cannot be treated as stochastically identical.

We next give an upper bound for $\eta_k.$ By virtue of (1.2),

\begin{eqnarray*}\mathbb{P}(X_{k+1}=0\,|\,X_k)\geq q^{\frac{n^2}{4}}\end{eqnarray*}

uniformly on $X_k,$ and hence

(2.7) \begin{eqnarray}\eta_k=\mathbb{P}(X_k>0)<\left(1-q^{\frac{n^2}{4}}\right)^k\qquad \forall\,k\in\mathbb{N}.\end{eqnarray}

The following result is a refinement of this naive estimate for the subcritical regime. By Markov’s inequality,

\begin{eqnarray*}\eta_k=\mathbb{P}(X_k\geq 1)\leq \mathbb{E}(X_k)=\mathbb{E}\Bigl(\sum_{x\in V_n}\mathcal{E}_k(x)\Bigr)=\sum_{x\in V_n}\xi_k(x).\end{eqnarray*}

This yields the following.

We remark that the result in the corollary will be further improved in Theorem 5.1 below, using a comparison to a branching process and known estimates for the extinction time of the latter. In particular, it turns out that if $np=1$ in (1.2), then

\begin{eqnarray*}\eta_k \leq 1-\Bigl(\frac{k}{k+2}\Bigr)^{i_0}\leq \frac{2i_0}{k+2}\qquad \forall\,k\in\mathbb{N},\end{eqnarray*}

and the bound is asymptotically tight (see the lower bound in Theorem 5.1 and also Corollary 5.1).

Proof of Proposition 2.3. We say that the distribution of $X_0$ is exchangeable if it is independent of a particular labeling of the network’s nodes; that is, $\mathbb{P}(B\subset A_0)=\mathbb{P}(\sigma B\subset A_0)$ for any permutation (one-to-one and onto relabeling) $\sigma:V_n\to V_n$ and any cluster of nodes $B\subset V_n.$

Fix any $x\in V_n.$ First, observe that if the distribution of $X_0$ is exchangeable, then

\begin{eqnarray*}\xi_k(x)=\mathbb{E}\bigl(\mathcal{E}_k(x)\bigr)=\frac{1}{n}\mathbb{E}\Bigl(\sum_{x\in V_n} \mathcal{E}_k(x)\Bigr)=\frac{1}{n}\mathbb{E}(X_k)\leq \frac{c^k}{n^2}\mathbb{E}(H_0)\end{eqnarray*}

by virtue of Corollary 2.1.

Consider now an auxiliary avalanche process where $A_0$ is chosen uniformly over subsets of $V_n$ of a given size m. We will denote the conditional law of the process $\mathbb{P}({\cdot}|A_0)$ by $\mathbb{P}_{m,x}$ when $x\in A_0$ and by $\mathbb{P}_{m,\overline x}$ otherwise. Then, in view of the result for exchangeable $X_0$ and Corollary 2.1,

\begin{eqnarray*}\frac{m(n-m)c^k}{n^2}&=&\frac{c^k}{n^2}\mathbb{E}(H_0)\geq \mathbb{P}\bigl(\mathcal{E}_k(x)=1\bigr)\\[2pt]&=&\frac{\binom{n-1}{m}}{\binom{n}{m}} \mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr)+\frac{\binom{n-1}{m-1}}{\binom{n}{m}} \mathbb{P}_{m,x}\bigl(\mathcal{E}_k(x)=1\bigr)\\[2pt] &\geq&\frac{\binom{n-1}{m}}{\binom{n}{m}} \mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr)=\frac{n-m}{n}\mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr),\end{eqnarray*}

which implies

(2.8) \begin{eqnarray}\mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr)\leq \frac{mc^k}{n}.\end{eqnarray}

Therefore,

(2.9) \begin{eqnarray}\nonumber\mathbb{P}_{j,x}\bigl(\mathcal{E}_{k+1}(x)=1\bigr)&=&\sum_{m=0}^{n-j}\binom{n-j}{m}\big(1-q^j\big)^mq^{j(n-j-m)}\mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr)\\\nonumber &\leq& c^k \sum_{m=0}^{n-j}\binom{n-j}{m}\big(1-q^j\big)^mq^{j(n-j-m)}\frac{m}{n}\\ &=&\frac{c^k}{n} \mathbb{E}_{j,x}(X_1)=\frac{c^k(n-j)\big(1-q^j\big)}{n} \leq \frac{c^{k+1}j(n-j)}{n^2},\end{eqnarray}

where $\mathbb{E}_{j,x}$ is the expectation with respect to the probability measure $\mathbb{P}_{j,x}.$ In particular, we have established that

(2.10) \begin{eqnarray}\mathbb{P}_{j,x}\bigl(\mathcal{E}_k(x)=1\bigr)\leq \frac{c^{k-1}(n-j)\big(1-q^j\big)}{n}\leq c^k\frac{j(n-j)}{n^2} \leq \frac{jc^k}{n}.\end{eqnarray}

Turning now to $\mathbb{P}_{j,\overline x}\bigl(\mathcal{E}_{k+1}(x)=1\bigr),$ write

\begin{eqnarray*}\mathbb{P}_{j,\overline x}\bigl(\mathcal{E}_{k+1}(x)=1\bigr)&=&\sum_{m=1}^{n-j-1}\binom{n-j-1}{m}\big(1-q^j\big)^mq^{j(n-j-m-1)}\mathbb{P}_{m,\overline x}\bigl(\mathcal{E}_k(x)=1\bigr)\\ && + \sum_{m=1}^{n-j-1}\binom{n-j-1}{m-1}\big(1-q^j\big)^{m-1}q^{j(n-j-m)}\mathbb{P}_{m,x}\bigl(\mathcal{E}_k(x)=1\bigr).\end{eqnarray*}

Using (2.10) along with (2.8), we obtain

(2.11) \begin{eqnarray}\nonumber\mathbb{P}_{j,\overline x}\bigl(\mathcal{E}_{k+1}(x)=1\bigr) &\leq&\frac{c^k}{n}\sum_{m=1}^{n-j-1}\binom{n-j-1}{m}\big(1-q^j\big)^mq^{j(n-j-m-1)}m\\\nonumber &&\qquad+ \frac{c^k}{n}\sum_{m=1}^{n-j-1}\binom{n-j-1}{m-1}\big(1-q^j\big)^{m-1}q^{j(n-j-m)}m\\ &=&\frac{c^k}{n} \mathbb{E}_{j,\overline x}(X_1)=\frac{c^k(n-j)\big(1-q^j\big)}{n}\leq \frac{c^{k+1}j(n-j)}{n^2}.\end{eqnarray}

The claim follows now from (2.9) and (2.11).

We next investigate the cumulative number of excited nodes in a subcritical regime. We will use the following lower bound for $1-q^i.$

Lemma 2.1. Under Condition 2.1,

\begin{equation*}1-q^i\geq \frac{di}{n}-\frac{c^2i^2}{2n^2}\end{equation*}

for all $i\in\mathbb{Z}_+.$

Proof. The lemma is trivial for $i=0,1.$ For $i\geq 2,$ using the Lagrange form of the second-order remainder in the Taylor series for $f(p)=(1-p)^i$ around zero,

\begin{eqnarray*}(1-p)^i=1-ip+\frac{i(i-1)}{2}p^2(1-p_*)^{i-2}\leq1-ip+\frac{i(i-1)}{2}p^2\leq 1- \frac{di}{n}+\frac{c^2i^2}{2n^2}\end{eqnarray*}

for some $p_*\in (0,p).$

Recall S from (1.3). We have the following.

Proposition 2.4. Suppose that Condition 2 holds with $c\in (0,1).$ Then

(2.12) \begin{eqnarray}\frac{\mathbb{E}(X_0)}{1-d}-\frac{3\mathbb{E}(X_0^2)}{n(1-c)^3}\leq \mathbb{E}(S) \leq \frac{\mathbb{E}(X_0)}{1-c}.\end{eqnarray}

Proof. To prove the proposition, we will first obtain a suitable lower bound for $\mathbb{E}(X_k).$ It follows from Lemma 2.1 that

(2.13) \begin{eqnarray}\nonumber\mathbb{E}(X_{k+1}|X_k)&\geq& (n-X_k)\Bigl(\frac{dX_k}{n}-\frac{c^2X_k^2}{2n^2}\Bigr)=dX_k-\frac{c^2X_k^2}{2n}-\frac{dX_k^2}{n}+\frac{c^2X_k^3}{2n^2}\\ \nonumber\\[-7pt] &\geq&dX_k-\frac{3cX_k^2}{2n},\end{eqnarray}

where we have used the fact that $d<c$ and $c^2<c.$ Using the identity in (2.4), we obtain

(2.14) \begin{eqnarray}\nonumber\mathbb{E}(X_k^2)&=&\mathbb{E}\bigl[(n-X_{k-1})\big(1-q^{X_{k-1}}\big)q^{X_{k-1}}+n^2\big(1-q^{X_{k-1}}\big)^2\bigr]\\\nonumber &\leq& \mathbb{E}\bigl[n(1-q^{X_{k-1}})+n^2\big(1-q^{X_{k-1}}\big)^2\bigr]\\ &\leq& \mathbb{E}\left[n\cdot \frac{c X_{k-1}}{n}+n^2 \frac{c^2 X_{k-1}^2}{n^2}\right]=\mathbb{E}\bigl(c X_{k-1}+c^2 X_{k-1}^2\bigr).\end{eqnarray}

Iterating,

\begin{eqnarray*}\mathbb{E}(X_k^2)&\leq&\mathbb{E}\bigl(c X_{k-1}+c^2 X_{k-1}^2\bigr)\leq \mathbb{E}\bigl(c X_{k-1}+c^3 X_{k-2} +c^4X_{k-2}^2\bigr)\\&\leq & c^{2k}\mathbb{E}(X_0^2)+\sum_{j=1}^k c^{2j-1}(X_{k-j})\leq c^{2k}\mathbb{E}(X_0^2)+\sum_{j=1}^k c^{k+j-1}\mathbb{E}(X_0)\\&\leq & c^{2k}\mathbb{E}(X_0^2)+\frac{c^k}{1-c}\mathbb{E}(X_0).\end{eqnarray*}

Therefore,

\begin{eqnarray*}E\mathbb{E}(X_{k+1})\geq d\mathbb{E}(X_k)-\frac{3c}{2n}\mathbb{E}\Bigl[c^{2k}\mathbb{E}(X_0^2)+\frac{c^k}{1-c}\mathbb{E}(X_0)\Bigr].\end{eqnarray*}

Iterating again, we obtain

\begin{eqnarray*}\mathbb{E}(X_{k+1})&\geq& d^{k+1}\mathbb{E}(X_0)-\frac{3c}{2n}\sum_{j=0}^k\Bigl[c^{2(k-j)}d^j\mathbb{E}(X_0^2)+\frac{c^{k-j}d^j}{1-c}\mathbb{E}(X_0) \Bigr]\\&\geq& d^{k+1}\mathbb{E}(X_0)-\frac{3c}{2n}\sum_{j=0}^k\Bigl[c^{2k+j}\mathbb{E}(X_0^2)+\frac{c^k}{1-c}\mathbb{E}(X_0) \Bigr]\\&\geq&d^{k+1}\mathbb{E}(X_0)-\frac{3c^{k+1}}{2n(1-c)}\bigl[c^k\mathbb{E}(X_0^2)+(k+1)\mathbb{E}(X_0) \bigr].\end{eqnarray*}

Using this bound along with the upper bound in Corollary 2.1, we obtain that

\begin{eqnarray*}d^k\mathbb{E}(X_0)-\frac{3c^k}{2n(1-c)}\bigl[c^k\mathbb{E}(X_0^2)+k\mathbb{E}(X_0) \bigr]\leq \mathbb{E}(X_k)\leq c^k\mathbb{E}(X_0).\end{eqnarray*}

Therefore, summing over all indexes k from zero to infinity,

\begin{eqnarray*}\frac{\mathbb{E}(X_0)}{1-d}-\frac{3}{2n(1-c)}\Bigl[\frac{1}{1-c^2}\mathbb{E}(X_0^2)+\frac{c}{(1-c)^2}\mathbb{E}(X_0) \Bigr]\leq \mathbb{E}\Bigl(\sum_{k=0}^\infty X_k\Bigr) \leq \frac{\mathbb{E}(X_0)}{1-c}.\end{eqnarray*}

Taking into account that $c< 1,$ $1-c^2<(1-c)^2,$ and $\mathbb{E}(X_0)\leq \mathbb{E}(X_0^2),$ we obtain the lower bound in the form given in the statement of the proposition.

We remark that though the constant

\begin{equation*}\frac{3\mathbb{E}(X_0^2)}{(1-c)^3}\end{equation*}

in front of $1/n$ in the correction term in the left-hand side of (2.12) is not optimal, the lower bound correctly captures the dependence of this term on n. The latter result is formally stated in Part (iii) of Theorem 3.1 below.

3. Poisson approximation and the size of the avalanche

The primary goal of this section is to study the asymptotic behavior of the cumulative size of the avalanche for a certain ensemble of comparable avalanche models. The underlying family of models is introduced in Equation (3.1) and Assumption 3.1 below, and the main result of this section is stated in Theorem 3.1. The secondary purpose of this section is to introduce a branching process approximation which will be used throughout the rest of the paper.

In the rest of the paper, along with a single network X, we will often consider a family of Markov chains $X^{(n)}=(X^{(n)}_k)_{k\in\mathbb{Z}_+},$ each governed by a transition kernel of the same type as in (1.2), namely

(3.1) \begin{eqnarray}\mathbb{P}_n(i,j)\,{:\!=}\,\mathbb{P}\bigl(X^{(n)}_{k+1}=j\bigl|X^{(n)}_k=i\bigr)=\binom{n-i}{j} (1-q_n^i)^j(q_n^i)^{n-i-j},\end{eqnarray}

for some $q_n\in (0,1)$ and all $i=0,1,\ldots,n,$ $j\leq n-i.$ In this definition we maintain the convention that $\binom{0}{0}=1$ and $\binom{0}{j}=0$ for all $i\in\mathbb{N}$ in (3.1). Therefore, all Markov chains in this collection eventually get absorbed at zero. Typically, we will impose the following comparability assumption on the family of avalanche models under consideration.

Some of our asymptotic estimates will be stated in terms of arbitrary numbers $c>0,$ $d>0,$ and $i_0\in \mathbb{N}$ satisfying the following condition. This condition is an analogue of Condition 2.1 for a family of networks that satisfies Assumption 3.1.

Condition 3.1. Let $p_n\in(0,1),$ $n\in\mathbb{N},$ be given, and consider a family of avalanche models $\{X^{(n)}:n\in\mathbb{N}\}$ with transition kernels defined in (3.1). Assume that Part (ii) of Assumption 3.1 is in force, and that furthermore there exist constants $c>0,$ $d\in (0,c),$ and $n_0\in\mathbb{N}$ such that

1. (i) $np_n \in [d,c]$ for all $n\geq n_0$;

2. (ii) if Part (i) of Assumption 3.1 holds and $\lambda>1,$ then $d>1$;

3. (iii) if Part (i) of Assumption 3.1 holds and $\lambda<1,$ then $c<1.$

It follows from (1.2) that under Assumption 3.1, for any $i\in\mathbb{N},$ $j\in\mathbb{Z}_+,$ and $k\in \mathbb{N},$

(3.2) \begin{eqnarray}\lim_{n\to\infty} \mathbb{P}\bigl(X^{(n)}_{k+1}=j\bigl|X^{(n)}_k=i\bigr)=e^{-\lambda i}\frac{(\lambda i)^j}{j!}.\end{eqnarray}

Let $Z^{(\lambda)}=\bigl(Z^{(\lambda)}_k\bigr)_{k\in \mathbb{Z}_+}$ be a Markov chain on $\mathbb{Z}_+$ with the sole absorption state at zero, a Poisson transition kernel

\begin{eqnarray*}\mathbb{P}\bigl(Z^{(\lambda)}_{k+1}=j\bigl|Z^{(\lambda)}_k=i\bigr)=e^{-\lambda i}\frac{(\lambda i)^j}{j!},\qquad i\in \mathbb{N},\ j\in\mathbb{Z}_+,\end{eqnarray*}

and the initial state $Z^{(\lambda)}_0=i_0,$ where $i_0$ is the common initial state of the Markov chains $X^{(n)}$ introduced in Assumption 3.1. We can assume without loss of generality that $Z^{(\lambda)}$ is a Galton–Watson branching process with a Poisson offspring distribution; that is,

(3.3) \begin{eqnarray}Z^{(\lambda)}_{k+1}=\sum_{j=1}^{Z^{(\lambda)}_k} Y^{(\lambda)}_{k,j}\end{eqnarray}

for a collection of independent Poisson random variables

\begin{eqnarray*}Y=\left\{Y^{(\mu)}_{k,j}:k\in \mathbb{Z}_+,j\in\mathbb{N},\mu>0\right\}\end{eqnarray*}

such that for all $i\in\mathbb{Z}_+,$

\begin{eqnarray*}\mathbb{P}\bigl(Y^{(\mu)}_{k,j}=i\bigr)=e^{-\mu}\frac{\mu^i}{i!}.\end{eqnarray*}

The sum in the right-hand side of (3.3) is assumed to be zero if $Z^{(\lambda)}_k=0$; that is, $Z^{(\lambda)}_k$ is formally defined for all $k\in\mathbb{Z}_+.$

The convergence in (3.2) implies the weak convergence of the sequence of Markov processes $X^{(n)}$ to the branching process $Z^{(\lambda)}$ as $n\to\infty$ [Reference Karr20]. To illustrate the functionality of the branching approximation, Figures 1 and 2 below provide plots of $\mathbb{E}(T)$ as a function of the initial state $i_0$ for $n=100$ and $n=1000,$ in each case for four values of the parameter $c=np$ concentrated around the theoretical phase transition value $c=1$ suggested by the branching approximation. Note that $\mathbb{E}(T)<\infty$ by virtue of (2.7).

Figure 1.

Plot of the function $f(i_0)=\mathbb{E}(T\,|\,X_0=i_0)$ for $n=100$ and several values of the parameter $c=np$ ranging from $c=0.9$ to $c=1.3.$

Figure 2.

Plot of the function $f(i_0)=\mathbb{E}(T\,|\,X_0=i_0)$ for $n=1000$ and several values of the parameter $c=np$ ranging from $c=0.9$ to $c=1.3.$

For an arbitrary $n\in\mathbb{N}$ and a given network $X^{(n)},$ let Q be $(n-1)\times (n-1)$ matrix with entries $Q(i,j)=\mathbb{P}\bigl(X^{(n)}_{k+1}=j\,\big|\,X^{(n)}_k=i\bigr).$ To simplify the notation, we do not indicate the index n for the avalanche duration T and the auxiliary matrix Q. To evaluate the expectation of T, we use the following standard Markov chain matrix calculation:

\begin{eqnarray*}\mathbb{E}(T)=\sum_{m=0}^\infty \mathbb{P}(T>m)=\sum_{m=0}^\infty Q^m e(i_0)=(I-Q)^{-1}e(i_0),\end{eqnarray*}

where $e\in\mathbb{R}^{n-1}$ is an $(n-1)$-vector with all entries equal to one and I is the $(n-1)$-dimensional unit matrix. To compute the inverse matrix in the above expression we used the packages ‘numpy’ and ‘decimal’ on Python 3.5 with the computation precision set to 400 decimal places.

An intuitive reason for the uniformly (on $i_0$) large values of $\mathbb{E}(T)$ and the high persistence of the avalanche in the supercritical regime, when n is large, is that the stochastic path of the Markov chain X is well approximated by a trajectory of a deterministic dynamical system that is locally Lipschitz, and consequently is quickly attracted to its unique (non-zero) global stable point (see Section 6 below for details). Heuristically, it appears that the Markov chain spends most of its time before the absorption being ‘trapped’ in a neighborhood of the stable point. The numerical simulations show that the phase transition in the avalanche model doesn’t occur at exactly $c=1$ for either $n=100$ or $n=1000.$ While the phase transition is fairly smooth for $n=100,$ it is considerably sharper and more conspicuous for $n=1000.$ Overall, one can conclude that the branching approximation gives a useful qualitative insight into the existence of an asymptotic phase transition in the avalanche model.

In what follows we will exploit the following explicit monotone coupling of $X^{(n)}$ with a branching process. For future convenience, we state the result in terms of a family of avalanche models rather than a single network. At the base of the construction is a standard coupling between a binomial B(n,p) and a Poisson$\bigl(-n\log(1-p)\bigr)$ random variable.

Proof of Proposition 3.1. We have

(3.4) \begin{eqnarray}-\log q_n=-\log(1-p_n)\leq -\log \Bigl(1-\frac{c}{n}\Bigr),\qquad \forall\ n\geq n_0.\end{eqnarray}

It is easy to check that for any $n\geq n_0$ and $i\in\mathbb{N},$

(3.5) \begin{eqnarray}-\log q_n\leq \frac{c}{n-1}\leq \frac{c}{n-i}.\end{eqnarray}

Indeed, if

\begin{equation*}f(x)=\frac{c}{x-1}+\log \Bigl(1-\frac{c}{x}\Bigr),\end{equation*}

then $\lim_{x\to +\infty}f(x)=0$ and for any $x>1,$

\begin{eqnarray*}f'(x)=-\frac{c}{(x-1)^2}+\frac{c}{x^2-cx}\leq -\frac{c}{(x-1)^2}+\frac{c}{x^2-x}=\frac{c}{x-1}\Bigl[-\frac{1}{x-1}+\frac{1}{x}\Bigr]<0,\end{eqnarray*}

implying that $f(x)>0$ for $x>1.$ Let

\begin{eqnarray*}\bigl\{Y^{(c,x,z)}_{n,k,j}:n\in\mathbb{N},k\in\mathbb{Z}_+,j\in\mathbb{N},x\in\mathbb{Z}_+,z\in\mathbb{Z}_+,x < n\bigr\}\end{eqnarray*}

be a collection of independent Poisson random variables such that

\begin{eqnarray*}\mathbb{P}\Big(Y^{(c,x,z)}_{n,k,j}=i\Big)=e^{-\frac{cz}{n-x}}\frac{\bigl(\frac{cz}{n-x}\bigr)^i}{i!}, \qquad i\in\mathbb{Z}_+.\end{eqnarray*}

Further, let $U=\bigl\{U^{(c,x)}_{n,k,j}:n\in\mathbb{N},k\in\mathbb{Z}_+,j\in\mathbb{N},x\in\mathbb{Z}_+,x< n\bigr\}$ be a collection of independent Bernoulli variables which is independent of the family of Poisson variables Y and satisfies

\begin{eqnarray*}\mathbb{P}\bigl(U^{(c,x)}_{n,k,j}=1\bigr)=\frac{1-q_n^x}{1-e^{-\frac{cx}{n-x}}}\qquad \textrm{and}\qquad \mathbb{P}\bigl(U^{(c,n)}_{n,k,j}=0\bigr)=\frac{q_n^x-e^{-\frac{cx}{n-x}}}{1-e^{-\frac{cx}{n-x}}}.\end{eqnarray*}

Finally, set

\begin{eqnarray*}B^{(c,x)}_{ n,k,j}=U^{(c,n)}_{n,k,j}\mbox 1_{\{{Y^{(c,x,x)}_{k,j}>0}\}}\end{eqnarray*}

and define a Markov chain of integer triples $(X^{(c,n)}_k,Q^{(c,n)}_k,Z^{(c,n)}_k)_{k\in\mathbb{Z}_+}$ through the initial condition $Z^{(c,n)}_0=Q^{(c,n)}_0=X^{(c,n)}_k=i_0$ and the recursion

(3.6) \begin{eqnarray}\left\{\begin{array}{ll}Z^{(c,n)}_{k+1}&=\sum_{j=1}^{n-x} \bigl(Y^{(c,x,x)}_{n,k,j}+Y^{(c,x,z-x)}_{n,k,j}\bigr),\\\mbox{}&\\Q^{(c,n)}_{k+1}&=\sum_{j=1}^{n-x} Y^{(c,x,x)}_{n,k,j},\\\mbox{}&\\X^{(c,n)}_{k+1}&=\sum_{j=1}^{n-x}B^{(c,x)}_{n,k,j},\end{array}\right.\qquad\textrm{if}\quad X^{(c,n)}_k=x,\ Z^{(c,n)}_k=z.\end{eqnarray}

By induction, $\mathbb{P}(X^{(c,n)}_k\leq Q^{(c,n)}_k\leq Z^{(c,n)}_k)=1$ for all $k\in\mathbb{N}.$ Furthermore, by our construction, $(X^{(c,n)}_k)_{k\in\mathbb{Z}_+}$ is distributed the same as $(X^{(n)}_k)_{k\in\mathbb{Z}_+}$ while $(Z^{(c,n)}_k)_{k\in\mathbb{Z}_+}$ is distributed the same as the branching process $(Z^{(c)}_k)_{k\in\mathbb{Z}_+}.$

Our next result concerns the cumulative size of the avalanche, namely the total number of excited sites created by the avalanche during its entire lifespan. Let

\begin{eqnarray*}S_n=\sum_{k=0}^\infty X^{(n)}_k.\end{eqnarray*}

Note that $\mathbb{P}(S_n<\infty)=1$ since $X^{(n)}$ is an irreducible Markov chain with a unique absorbing state at zero. The following theorem complements the bounds provided by Proposition 2.4 for a single network with a fixed $n\in\mathbb{N}.$ The theorem relates asymptotic characteristics of $X^{(n)}$ to their counterparts for the limiting branching process $Z^{(\lambda)}.$

For $\mu>0,$ let $S^{(\mu)}=\sum_{k=0}^\infty Z^{(\mu)}_k.$ By the Otter–Dwass theorem, the limiting distribution in (3.7) is the distribution of $S^{(\lambda)}$ [Reference Dwass17]. Similarly to other models using an approximation by the Poisson branching process, the distribution tails of the cumulative size of the underlying population (avalanche in our case) at the critical regime obey a power law. Indeed, (3.7) and Stirling’s formula imply that when $\lambda=1,$ for large values of j and n, $\mathbb{P}(S_n=j)$ is well approximated by

\begin{equation*}\frac{i_0}{j}\frac{(\lambda j)^{j-i_o}}{(j-i_0)!}e^{-\lambda j}\sim i_0\sqrt{\frac{1}{2\pi}}j^{-3/2}.\end{equation*}

Proof of Theorem 3.1.

1. (i) If $\lambda<1$, the result in (i) follows from Proposition 2.4. If $\lambda\geq 1,$ then a version of Fatou’s lemma for weakly convergent sequences implies that for any $J\in\mathbb{N},$

   \begin{eqnarray*}\liminf_{n\to\infty} \mathbb{E}\left(\sum_{k=0}^J X^{(n)}_k\right)\geq \mathbb{E}\left(\sum_{k=0}^J Z^{(\lambda)}_k\right)=i_0\sum_{k=0}^J \lambda^k,\end{eqnarray*}
   and the result follows by taking J to infinity. We remark in passing that estimates similar to (2.12) show that in fact
   \begin{equation*}\lim_{n\to\infty} \mathbb{E}\left(\sum_{k=0}^J X^{(n)}_k\right)=i_0\sum_{k=0}^J \lambda^k.\end{equation*}

2. (ii) Let $c>0$ and $d>0$ be as in Condition 3.1. Assume first that $\lambda\in(0,1).$ A simple argument can be given to prove the result in this case. To prove the convergence of $S_n$ to $S^{(\lambda)}$ we will consider exponential generating functions $\mathbb{E}(e^{-\alpha S_n})$ and $\mathbb{E}(e^{-\alpha S^{(c)}}),$ $\alpha>0,$ $c>0,$ and use the inequality $e^{-\alpha x}-e^{-\alpha y}\leq \alpha(y-x)$, which is true for any $y>0$ and $x\in(0,y).$ It follows from (3.6) that

   \begin{eqnarray*}0\leq \mathbb{E}\big(e^{-\alpha S_n}\big)-\mathbb{E}\big(e^{-\alpha S^{(c)}}\big)\leq \alpha \mathbb{E}\big(S^{(c)}-S_n\big).\end{eqnarray*}
   By Proposition 2.4, for any $n\geq n_0,$
   \begin{eqnarray*}0\leq \limsup_{n \to\infty}\bigl[\mathbb{E}(e^{-\alpha S_n})-\mathbb{E}(e^{-\alpha S^{(c)}})\bigr]\leq \alpha\Bigl(\frac{i_0}{1-c}-\frac{i_0}{1-d}\Bigr),\end{eqnarray*}
   which yields the result since the parameters c and d can be chosen arbitrarily close to $\lambda.$

   Assume now that $\lambda=1.$ Without loss of generality we can assume that $c>1.$ Let $A_c=\bigl\{\lim_{k\to\infty} Z^{(c)}_k=0\bigr\}$ be the event of extinction for the branching process $\bigl(Z^{(c)}_k\bigr)_{k\in\mathbb{Z}_+}.$ It follows from Proposition 3.1 that for any integer $m\geq i_0,$

   (3.8) \begin{eqnarray}\mathbb{P}(S_n>m)\leq \mathbb{P}\bigl(S^{(c)}>m;A_c\bigr)+\mathbb{P}\bigl(\overline A_c\bigr)=\sum_{j>m}\frac{i_0}{j}\frac{(c j)^{j-i_o}}{(j-i_0)!}e^{-c j}+\mathbb{P}\bigl(\overline A_c\bigr),\end{eqnarray}
   where $\overline A_c$ is the complement event $A_c=\bigl\{\lim_{k\to\infty} Z^{(c)}_k=+\infty\bigr\},$ and the second identity is an instance of the Otter–Dwass theorem for supercritical branching processes (see Theorem 1 in [Reference Dwass17]). Recall that $\mathbb{P}(A_c)$ is the smallest nonnegative root of the equation $e^{t(c-1)}=t$ [Reference Athreya and Ney5, p. 7], and hence $\mathbb{P}(\overline A_c)$ is a continuous function of the parameter c. Therefore, by first letting n go to infinity and then letting c approach $\lambda=1,$ we obtain
   (3.9) \begin{eqnarray}\limsup_{n\to\infty}\mathbb{P}(S_n>m)\leq \sum_{j>m}\frac{i_0}{j}\frac{(j)^{j-i_o}}{(j-i_0)!}e^{-j}.\end{eqnarray}
   On the other hand, Fatou’s lemma implies that for any $J\in\mathbb{N},$
   (3.10) \begin{align}\liminf_{n\to\infty} \mathbb{P}(S_n>m)\geq \liminf_{n\to\infty} \mathbb{P}\Bigl(\sum_{k=0}^J X^{(n)}_k>m\Bigr)\geq\mathbb{P}\Bigl(\sum_{k=0}^J Z^{(1)}_k>m\Bigr).\end{align}
   Letting first n and then J go to infinity, we obtain that
   (3.11) \begin{eqnarray}\liminf_{n\to\infty}\mathbb{P}(S_n>m)\geq \mathbb{P}\bigl(S^{(1)}>m\bigr)=\sum_{j>m}\frac{i_0}{j}\frac{(j)^{j-i_o}}{(j-i_0)!}e^{-j}.\end{eqnarray}
   Combining this estimate with (3.9) completes the proof of Part (ii) for $\lambda=1.$

3. (iii) The proof is similar to that of Part (ii) for $\lambda=1.$ More precisely, (3.8) and (3.10) with $Z^{(1)}$ replaced by $Z^{(\lambda)}$ remain correct for any $m>i_0,$ $c>\lambda,$ and $J\in\mathbb{N}.$ By first letting n go to infinity and then letting c approach $\lambda$ in (3.8), we obtain the following counterpart of (3.9):

   (3.12) \begin{eqnarray}\limsup_{n\to\infty}\mathbb{P}(S_n>m)\leq \sum_{j>m}\frac{i_0}{j}\frac{(j\lambda)^{j-i_o}}{(j-i_0)!}e^{-j\lambda}+1-\alpha_\lambda.\end{eqnarray}
   By letting first n and then J go to infinity in (3.11), we obtain from the Otter–Dwass theorem that
   \begin{eqnarray*}\liminf_{n\to\infty}\mathbb{P}(S_n>m)\geq \mathbb{P}\bigl(S^{(\lambda)}>m\bigr)=\sum_{j>m}\frac{i_0}{j}\frac{(j\lambda)^{j-i_o}}{(j-i_0)!}e^{-j\lambda}+1-\alpha_\lambda.\end{eqnarray*}
   Combining this estimate with (3.12) completes the proof of Part (iii) of the theorem.

4. (iv) Using the Lagrange form of the second-order remainder in the Taylor series for the function $f(p)=(1-p)^i$ around zero, we obtain that for all $n\in\mathbb{N}$ and $i\in\mathbb{N},$

   \begin{eqnarray*}q_n^i=(1-p_n)^i=1-ip_n+\frac{i(i-1)}{2}p_n^2(1-\beta_{n,i})^{i-2}\end{eqnarray*}
   for some $\beta_{n,i}\in (0,p_n).$ Therefore,
   (3.13) \begin{eqnarray}\nonumber\mathbb{E}(X^{(n)}_{k+1})&=& \mathbb{E}\left[\big(n-X^{(n)}_k\big)\left(1-q_n^{X^{(n)}_k}\right)\right]\\ &=&np_n\mathbb{E}\big(X^{(n)}_k\big) +\frac{np_n^2}{2}\mathbb{E}\left[X^{(n)}_k\big(X^{(n)}_k-1\big)\right]\end{eqnarray}
   \begin{eqnarray}\nonumber&&-p_n\mathbb{E}[\big(X^{(n)}_k\big)^2]+\frac{p_n^2}{2}\mathbb{E}\left[\big(X^{(n)}_k\big)^2\big(X^{(n)}_k-1\big)\bigg(1-\beta_{n,X^{(n)}_k}\bigg)^{X^{(n)}_k-2}\right].\end{eqnarray}
   It follows from the coupling construction given by Proposition 3.1 that
   \begin{eqnarray*}\sum_{k=0}^\infty \mathbb{E}\bigl[(X^{(n)}_k)^3\bigr]<\sum_{k=0}^\infty \mathbb{E}\bigl[(Z^{(c)}_k)^3\bigr]<\infty\end{eqnarray*}
   for any $c\in (\lambda,1)$ and $n\geq n_0.$ Hence, summing both sides of (3.13) from $k=1$ to infinity, we obtain
   \begin{eqnarray*}\mathbb{E}(S_n)&=&\frac{i_0}{1-np_n}+\frac{np_n^2}{2(1-np_n)}\mathbb{E}\left[\sum_{k=0}^\infty X^{(n)}_k\big(X^{(n)}_k-1\big)\right]\\ &&\qquad-\frac{p_n}{1-np_n}\mathbb{E}\left[\sum_{k=0}^\infty \big(X^{(n)}_k\big)^2\right]+o(1/n).\end{eqnarray*}
   Therefore, by the dominated convergence theorem (here we use again Proposition 3.1, which shows that $X^{(n)}_k$ is stochastically dominated by $Z^{(c)}_k$),
   (3.14) \begin{eqnarray}\nonumber&&\lim_{n\to\infty}n\Bigl[\frac{i_0}{1-c}-\mathbb{E}(S_n)\Bigr]\\\nonumber &&\qquad=-\frac{\lambda^2}{2(1-\lambda)}\mathbb{E}\left[\sum_{k=0}^\infty \big(Z^{(\lambda)}_k\big)^2-\frac{i_0}{1-\lambda}\right]+\frac{\lambda}{1-\lambda}\mathbb{E}\left[\sum_{k=0}^\infty \big(Z^{(\lambda)}_k\big)^2\right]\\ &&\qquad=\frac{i_0\lambda^2}{2(1-\lambda)^2}+\frac{2\lambda-\lambda^2}{2(1-\lambda)}\mathbb{E}\left[\sum_{k=0}^\infty \big(Z^{(\lambda)}_k\big)^2\right].\end{eqnarray}
   Furthermore,
   \begin{eqnarray*}\mathbb{E}\left[\sum_{k=0}^\infty \big(Z^{(\lambda)}_k\big)^2\right]=i_0^2+\mathbb{E}\left[\sum_{k=1}^\infty \left[\lambda Z^{(\lambda)}_{k-1}+ \lambda^2 \big(Z^{(\lambda)}_{k-1}\big)^2\right]\right],\end{eqnarray*}
   which implies the following (the sum is finite, for instance, because it is dominated by a finite second moment of the Borel–Tanner distribution of $S^{(\lambda)}$):
   \begin{eqnarray*}\mathbb{E}\left[\sum_{k=0}^\infty \big(Z^{(\lambda)}_k\big)^2\right]=\frac{\lambda i_0}{(1-\lambda)(1-\lambda^2)}+\frac{i_0^2}{1-\lambda^2}.\end{eqnarray*}
   Substituting this identity into (3.14) yields the result in Part (iv) of the theorem.

4. Spread to a non-zero fraction of the network

In this section we are concerned with the question of whether an avalanche initiated by just a few excited nodes has substantial potential to spread to a large fraction of the network. The results for the supercritical regime are stated in Theorems 4.1 and 4.2, while the critical and subcritical regimes are addressed in Theorem 4.4.

First we will consider a single network with given parameters n and p. For an arbitrary real number $J>0$ we define

(4.1) \begin{eqnarray}h_J(i)=\mathbb{P}\Bigl(\max_{k\in\mathbb{Z}_+} X_k \geq J\Bigl|X_0=i\Bigr)=\mathbb{P}\bigl(X_{T_J} \geq J|X_0=i\bigr),\end{eqnarray}

where

(4.2) \begin{eqnarray}T_J=\min\bigl\{k\in\mathbb{N}:X_k=0\ \mbox{or}\ X_k \geq J\bigr\}.\end{eqnarray}

We begin with the supercritical regime, namely the case when $np>1.$ Consequently, without loss of generality we can assume that $d>1$ in Condition 2.1. Let $\varepsilon_d>0$ be a positive constant such that

(4.3) \begin{eqnarray}\frac{1-e^{-d\varepsilon}}{\varepsilon}(1-\varepsilon)>1 \qquad \forall\ \varepsilon\in (0,\varepsilon_d].\end{eqnarray}

Note that such $\varepsilon_d$ exists because

\begin{equation*}\lim_{\varepsilon\to 0} \frac{1-e^{-d\varepsilon}}{\varepsilon}(1-\varepsilon)=d>1.\end{equation*}

Further, for $\mu\in (0,\infty)\backslash\{1\},$ let $\alpha_\mu\neq 1$ denote the unique solution in $(0,\infty)\backslash\{1\}$ of the fixed point equation

(4.4) \begin{eqnarray}\alpha_\mu=e^{-(1-\alpha_\mu)\mu},\qquad \alpha_\mu\neq 1.\end{eqnarray}

We remark that $\alpha_\mu<\mu\alpha_\mu<1$ if $\mu>1,$ and $\alpha_\mu>\mu\alpha_\mu>1$ if $\mu\in (0,1).$ This is true because (4.4) is equivalent to $\mu\alpha_\mu e^{-\mu\alpha_\mu}=\mu e^{-\mu},$ and the function $f(x)=xe^{-x}$ has a unique local maximum at $x=1.$

Observe that the right-hand side of (4.4) is

\begin{equation*}\mathbb{E}\Bigl(\alpha_\mu^{Z^{(\mu)}_{k+1}}\Bigl|Z^{(\mu)}_k=1\Bigl),\end{equation*}

and hence the sequence

\begin{equation*}\Bigl(\alpha_\mu^{Z^{(\mu)}_k}\Bigr)_{k\in\mathbb{Z}_+}\end{equation*}

is a martingale with respect to its natural filtration. If $\mu>1,$ then $\alpha_\mu$ is the extinction probability of the supercritical Poisson branching process $Z^{(\mu)}.$ On the other hand, if $\mu<1$, then $1/\alpha_\mu$ is the extinction probability of the dual supercritical process $Z^{(\mu\alpha_\mu)}.$

The main technical result of this section is the following proposition.

Proof.

1. (a) First, we choose a real constant $\gamma>0$ in such a way that

   \begin{eqnarray*}q^i\leq \Bigl(1-\frac{d}{n}\Bigr)^i\leq 1-\frac{\gamma i}{n}\end{eqnarray*}
   for any integer $i\in [1,n\varepsilon_d).$ Since
   \begin{eqnarray*}q^i\leq \Bigl(1-\frac{d}{n}\Bigr)^i=\Bigl[\Bigl(1-\frac{d}{n}\Bigr)^n\Bigr]^{\frac{i}{n}}\leq e^{-\frac{di}{n}},\end{eqnarray*}
   it suffices to find $\gamma >0$ such that $e^{-dx}\leq 1-\gamma x$ for any $x\in [0,\varepsilon_d).$ To this end, for a fixed $\gamma>0$ let $f_\gamma(x)=1-\gamma x-e^{-dx}.$ Then $f_\gamma(0)=0,$ $f_\gamma^{\prime}(x)=-\gamma+de^{-dx},$ and $f_\gamma^{\prime\prime}(x)=-d^2e^{-dx}<0.$ Thus $f_\gamma(x)>0$ for any $x\in (0,\varepsilon_d)$ provided that
   \begin{eqnarray*}f_\gamma^{\prime}(0)=-\gamma+d>0\qquad \textrm{and}\qquad f_\gamma(\varepsilon_d)=1-\gamma\varepsilon_d -e^{-d\varepsilon_d}\geq 0.\end{eqnarray*}
   Since
   \begin{equation*}\frac{1-e^{-d\varepsilon_d}}{\varepsilon_d}<d,\end{equation*}
   we can put
   (4.6) \begin{eqnarray}\gamma=\frac{1-e^{-d\varepsilon_d}}{\varepsilon_d}.\end{eqnarray}
   Then, for any $\rho\in (0,1)$ and integer $i\in [1,n\varepsilon_d),$
   (4.7) \begin{eqnarray}\nonumber\mathbb{E}\bigl(\rho^{X_{k+1}}\bigr|X_k=i\bigr)&=& \bigl(\rho (1-q^i)+1\cdot q^i\bigr)^{n-i}=\bigl(\rho +(1-\rho)\cdot q^i\bigr)^{n-i}\\\nonumber&\leq&\Bigl(\rho +(1-\rho)\cdot \Bigl(1-\frac{\gamma i}{n}\Bigr)\Bigr)^{n-i}=\Bigl(1 -(1-\rho)\frac{\gamma i}{n}\Bigr)^{n-i}\\&\leq& e^{-(1-\rho)\gamma(1-\varepsilon_d) i}.\end{eqnarray}
   Thus, we can set $\rho$ to be the unique solution of the fixed point equation
   (4.8) \begin{eqnarray}\rho= e^{-(1-\rho)\gamma (1-\varepsilon_d)}.\end{eqnarray}
   Note that $\rho=\alpha_\mu$ with
   \begin{equation*}\mu=\gamma (1-\varepsilon_d)=\frac{1-e^{-d\varepsilon_d}}{\varepsilon_d} (1-\varepsilon_d).\end{equation*}
   The limit result in (4.5) follows immediately from (4.8) and (4.6).

2. (b) By Proposition 3.1 the process $Z^{(c)}$ stochastically dominates X. Therefore, taking into account that $\alpha_c<1.$ we obtain

   \begin{eqnarray*}\mathbb{E}\bigl(\alpha_c^{X_{k+1}}|X_k=i\bigr)\geq \mathbb{E}\left(\alpha_c^{Z^{(c)}_{k+1}}|Z^{(c)}_k=i\right)=\alpha_c^i.\end{eqnarray*}
   The proof of the proposition is complete.

Doob’s optional stopping theorem implies that for any $\varepsilon\in (0,\varepsilon_d)$ and $i\in [1,n\varepsilon),$

\begin{eqnarray*}\rho^i\geq \mathbb{E}\bigl(\rho^{X_{T_{n\varepsilon}}}|X_0=i\bigr)\geq h_{n\varepsilon}(i)\rho^{n\varepsilon}+\bigl(1-h_{n\varepsilon}(i)\bigr)\end{eqnarray*}

and

\begin{eqnarray*}\alpha_c^i\leq \mathbb{E}\bigl(\alpha_c^{X_{T_{n\varepsilon}}}|X_0=i\bigr)\leq h_{n\varepsilon} (i)\alpha_c^{n\varepsilon}+\bigl(1-h_{n\varepsilon} (i)\bigr).\end{eqnarray*}

This yields the following result.

Theorem 4.1. Suppose that Condition 2.1 is satisfied with $d>1.$ Then

\begin{eqnarray*}\frac{1-\rho^i}{1-\rho^{n\varepsilon}}\leq h_{n\varepsilon}(i)\leq \frac{1-\alpha_c^i}{1-\alpha_c^{n \varepsilon}}\end{eqnarray*}

for all $\varepsilon\in (0,\varepsilon_d)$ and $i\in [1, n \varepsilon).$

We next consider the asymptotic behavior of the avalanche model when the network size approaches infinity. Similarly to (4.1) and (4.2), for any $J>0$ and $n\in\mathbb{N}$ we define

(4.9) \begin{eqnarray}h^{(n)}_J(i)=\mathbb{P}\Bigl(\max_{k\in\mathbb{Z}_+} X^{(n)}_k \geq J\Bigl|X^{(n)}_0=i\Bigr)=\mathbb{P}\bigl(X^{(n)}_{T^{(n)}_J} \geq J\bigl|X^{(n)}_0=i\bigr),\end{eqnarray}

where

(4.10) \begin{eqnarray}T^{(n)}_J=\min\bigl\{k\in\mathbb{N}:X^{(n)}_k=0\ \mbox{or}\ X^{(n)}_k \geq J \bigr\}.\end{eqnarray}

Suppose that Assumption 3.1 holds (and consequently Condition 3.1 holds) with $\lambda>1.$ Theorem 4.1 then implies that for any constants $i\in \mathbb{N},$ $\varepsilon\in (0,\varepsilon_d),$ and a function $\psi:\mathbb{N}\to\mathbb{N}$ such that

(4.11) \begin{eqnarray}\lim_{n\to\infty} \psi(n)=+\infty\qquad \textrm{and}\qquad \limsup_{n\to\infty} \frac{\psi(n)}{n}=0,\end{eqnarray}

we have

\begin{equation*}1-\rho^i \leq \liminf_{n\to\infty} h^{(n)}_{\psi(n)}(i)\leq \limsup_{n\to\infty} h^{(n)}_{\psi(n)}(i)\leq 1-\alpha_c^i,\end{equation*}

where $\rho=\rho(\varepsilon),$ as defined in the statement of Part (a) of Proposition 4.1. Because of the second condition in (4.11), we can choose the constant $\varepsilon>0$ to be as small as we wish. Therefore, by virtue of (4.5),

\begin{eqnarray*}1-\alpha_d^i \leq \liminf_{n\to\infty} h^{(n)}_{\psi(n)}(i)\leq \limsup_{n\to\infty} h^{(n)}_{\psi(n)}(i)\leq 1-\alpha_c^i.\end{eqnarray*}

Since the constants c and d can be chosen arbitrarily close to $\lambda,$ we arrive at the following theorem.

Theorem 4.2. Suppose that Assumption 3.1 is satisfied with $\lambda>1.$ Then, for any integer $i\in \mathbb{N}$ and a function $\psi:\mathbb{N}\to \mathbb{N}$ which satisfies (4.11), we have

(4.12) \begin{eqnarray}\lim_{n\to\infty} h^{(n)}_{\psi(n)}(i)=1-\alpha_\lambda^i.\end{eqnarray}

We remark that a similar result for a frequency-dependent Wright–Fisher model has been obtained in [Reference Chumley, Aydogmus, Matzavinos and Roitershtein13]. An anonymous reviewer pointed out that a limit similar to (4.12) should exist also for $h^{(n)}_{n\varepsilon}(i)$ with arbitrary, or perhaps at least large enough, $\varepsilon\in (0,1).$ We state this as a conjecture.

We now turn to the study of the maximal number of excited sites in the subcritical and critical regimes.

Remark 4.1. Let $M^{(\lambda)}=\max_{k\in\mathbb{Z}_+} Z^{(\lambda)}_k$ and $M_n=\max_{k\in\mathbb{Z}_+} X^{(n)}_k.$ The estimates on the right-hand side of (4.13) and (4.14) are classical upper bounds for the probability $\mathbb{P}\bigl(M^{(\lambda)}\geq m\bigl| Z^{(\lambda)}_0=i\bigr)$ [Reference Lindvall27]. Note that the estimates are not trivial, in the sense that, in general, $X^{(n)}$ is not dominated by the limiting branching process $Z^{(\lambda)}$, because it is possible that $np_n>\lambda.$ However, since both $M^{(\lambda)}$ and $M_n$ are a priori finite with probability one, (3.2) implies that $M_n$ converges to $M^{(\lambda)}$ in distribution as $n\to\infty.$ Hence one can expect that the bounds are meaningful for the avalanche model when n is large. The bound in (4.14) is known to be asymptotically accurate as $m\to\infty;$ that is,

\begin{equation*}\lim_{m\to\infty} \mathbb{P}\bigl(M^{(1)}\geq m\bigl| Z^{(1)}_0=i\bigr)=\frac{i}{m}\end{equation*}

(see [Reference Lindvall27]). For the subcritical process, Theorem 2 in [Reference Nerman32] suggests that the correct order of $\mathbb{P}\bigl(M^{(\lambda)}\geq m\bigl| Z^{(\lambda)}_0=i\bigr)$ as $m\to\infty$ is $m^{-1}\alpha_\lambda^{-m},$ up to a constant that depends on i. Before we prove Theorem 4.4 we state a direct consequence for our model of some well-known results for branching processes in the critical and subcritical regimes. The first part is an implication of a result in [Reference Lindvall27] mentioned in Remark 4.1, the second part can be derived from a result of [Reference Athreya3], and the third one follows from Theorem 2 in [Reference Nerman32]; all three are based on a comparison to a branching process invoking Proposition 3.1. We continue to use the notation introduced in the above remark for maxima of a branching process.

We remark that we chose the sequence $\varepsilon_m$ in the statement of the proposition to be the same in all three cases exclusively for notational convenience.

We now return to Theorem 4.4.

Proof of Theorem 4.4.

1. (a) Let $n_0\in\mathbb{N}$ and $d\in (0,\lambda)$ satisfy Condition 3.1. Fix any $\varepsilon>0$ and, similarly to (4.6), let

   (4.15) \begin{eqnarray}\gamma=\gamma(\varepsilon)=\frac{1-e^{-d\varepsilon}}{\varepsilon}.\end{eqnarray}
   Notice that $\gamma<d<1.$ Similarly to (4.7), for any $n\geq n_0,$ $\rho>1$, and integer $i\in [1,n\varepsilon),$
   \begin{eqnarray*}\mathbb{E}\bigl(\rho^{X^{(n)}_{k+1}}\bigr|X^{(n)}_k=i\bigr)\leq\Bigl(1 +(\rho-1)\frac{\gamma i}{n}\Bigr)^{n-i}\leq e^{(\rho-1)\gamma i}.\end{eqnarray*}
   Recall (4.4) and set $\rho=\alpha_\gamma.$ Thus
   \begin{equation*}\mathbb{E}\left(\rho^{X^{(n)}_{k+1}}\biggr|X^{(n)}_k=i\right)\leq \rho^i\end{equation*}
   for any $n\geq n_0$ and $\in [1,n\varepsilon).$ Doob’s optional theorem implies that
   \begin{eqnarray*}\rho^i\geq \mathbb{E}\left(\rho^{X^{(n)}_{T^{(n)}_m}}\Biggl|X^{(n)}_0=i\right)\geq h^{(n)}_m(i)\rho^m+\bigl(1-h^{(n)}_m(i)\bigr).\end{eqnarray*}
   Therefore,
   \begin{eqnarray*}h^{(n)}_m(i)\leq \frac{\rho^i-1}{\rho^m-1}=\frac{\alpha_\gamma^i-1}{\alpha_\gamma^m-1}\qquad \forall\ n\geq n_0.\end{eqnarray*}
   By taking $\varepsilon$ to zero in (4.15) we obtain
   (4.16) \begin{eqnarray}\limsup_{n\to\infty} h^{(n)}_m(i)\leq \frac{\alpha_d^i-1}{\alpha_d^m-1}.\end{eqnarray}
   Since in this argument d is an arbitrary number in $(0,\lambda),$ the result of Part (a) follows from taking d to $\lambda$ in the above inequality.

2. (b) Observe that (4.16) continues to hold for $\lambda=1$ and any $d\in (0,1).$ Furthermore, $\lim_{d\to 1} \alpha_d=1,$ and hence,

   \begin{equation*}\limsup_{n\to\infty} h^{(n)}_m(i)\leq \lim_{\alpha\to 1}\frac{\alpha^i-1}{\alpha^m-1}=\frac{i}{m}\end{equation*}
   by L’Hospital’s rule. The proof of the proposition is complete.

5. Duration of the avalanche

The goal of this section is to evaluate distribution tails of the avalanche’s duration. The expected value of the duration is discussed in some detail elsewhere [Reference Longini28, Reference Rastegar and Roitershtein35] using a mixture of numerical and computational approaches. The main result here is Theorem 5.1, which provides estimates for a single network. Some consequences for a network ensemble satisfying Assumption 3.1 are drawn in Corollaries 5.1 and 5.2 at the end of the section. The basic idea of the proofs is to compare the avalanche model to a branching process during a time frame which on one hand is large enough for an asymptotic pattern to emerge, and on the other hand is sufficiently small so that with an asymptotically overwhelming probability, the paths of the avalanche Markov chain X and a coupled branching process will not diverge until the end. The proofs rely on asymptotically tight estimates for the branching process given in [Reference Agresti1].

For $c\in (0,1)$ let

(5.1) \begin{eqnarray}s(s)=\frac{2-c}{c} \qquad \mbox{and} \qquad r(c)=\frac{ce^{-c}}{e^{-c}-(1-c)}.\end{eqnarray}

Recall T from (1.3) and $\alpha_c$ from (4.4).

Note that the right-hand side in (5.2) can be optimized with respect to J, and the optimal parameter value $J_{opt}$ is a solution of the equation

\begin{equation*}\frac{3cmJ}{n}-\frac{(\alpha_c^{i_0}-1)\ln \alpha_c}{(\alpha_c^J-1)^2}=0,\end{equation*}

implying that $J_{opt}\cong const\cdot (n/m)^3$ for large n.

The above bounds for the distribution function of T originate from their counterparts for the limiting branching process; see Lemma 5.2 below. The latter estimates are borrowed from [Reference Agresti1]. We remark that in the supercritical regime $c>1,$ a similar result for a frequency-dependent Wright–Fisher model has been proved in [Reference Chumley, Aydogmus, Matzavinos and Roitershtein13] (see Theorem 3.9 there). By taking n to infinity in the conclusions of the theorem, one can obtain tight asymptotic bounds for the avalanche model under Assumption 3.1; see Corollaries 5.1 and 5.2 below for details.

Proof of Theorem 5.1. Let $d_{TV}(X,Y)$ denote the total variation distance between the distributions of the random variables X and Y. That is, if X and Y are both non-negative and integer-valued, $d_{TV}(X,Y)=\frac{1}{2}\sum_{n=0}^\infty |\mathbb{P}(X=n)-\mathbb{P}(Y=n)|.$ By the coupling inequality, $d_{TV}(X,Y)\leq \mathbb{P}(X\neq Y).$ Furthermore, there exists a maximal coupling, that is, a random pair $\bigl(\widetilde X,\widetilde Y\bigr)$ such that $\widetilde X$ is distributed the same as X, $\widetilde Y$ is distributed the same as Y, and $\mathbb{P}\bigl(\widetilde X\neq \widetilde Y)=d_{TV}(X,Y)$ [Reference Thorisson37].

We will use the following inequalities. For the first claim, see, for instance, Theorem 4 and the subsequent Remark 1.1.4 in [Reference Chatterjee, Diaconis and Meckes12], and for the second one, see Theorem 1.C(i) in [Reference Barbour, Holst and Janson9].

Using the above results, we can construct a coupling of the Markov chain $(X_k)_{k\in \mathbb{Z}_+}$ and the limiting branching process $(Z_k)_{k\in \mathbb{Z}_+}$ as follows. The resulting process $\bigl(\widetilde X_k,\widetilde Z_k\bigr)_{k\in \mathbb{Z}_+}$ is a Markov chain. Suppose that the random pairs $\bigl(\widetilde X_t,\widetilde Z_t\bigl)$ have been sampled for all $t\leq k$ and that $\widetilde X_t=\widetilde Z_t$ for all $t\leq k.$ Let i be the common value of $\widetilde Z_k$ and $\widetilde X_k.$ Sample the next pair $\bigl(\widetilde X_{k+1},\widetilde Z_{k+1}\bigr)$ using the maximal coupling for $X_{k+1}$ under the conditional distribution $P\mathbb{P}(X_{k+1}\in \cdot\,|X_k=i)$ and $Z^{(c)}_{k+1}$ under the conditional distribution $\mathbb{P}\bigl(Z^{(c)}_{k+1}\in \cdot\,\bigl|Z^{(c)}_k=i\bigr).$ After the random time $\tau\,{:\!=}\,\min\bigl\{k\in\mathbb{N}:\widetilde X_k\neq \widetilde Z_k\bigr\},$ sample $\bigl(\widetilde Z_t\bigr)_{t\geq \tau}$ and $\bigl(\widetilde X_t)_{t\geq \tau}$ independently. Using Lemma 5.1 and at first approximating $\widetilde X_{k+1}$ by a Poisson random variable with parameter $(n-i)(1-q^i),$ we obtain that for $c>1,$

\begin{align*}\nonumber\mathbb{P}\bigl(\widetilde X_{k+1}\neq \widetilde Z_{k+1}\bigl|\widetilde X_k=\widetilde Z_k=i\bigr)&\leq \frac{1}{2}\big(1-q^i\big)+\frac{1}{\sqrt{c i}}\left|(n-i)\big(1-q^i\big)-c i\right|\\&\leq \frac{ci}{2n}+\frac{1}{\sqrt{c i}}\bigl[ci-(n-i)(1-q^i)\bigr].\end{align*}

Using the bound in Lemma 2.1 we conclude that in the case when $c>1,$

(5.3) \begin{eqnarray}\mathbb{P}\bigl(\widetilde X_{k+1}\neq \widetilde Z_{k+1}\bigl|\widetilde X_k=\widetilde Z_k=i\bigr)\leq \frac{ci}{2n}+\frac{c^{1/2}i^{3/2}}{2n}+\frac{c^{3/2}i^{3/2}}{2n}\leq \frac{3c^{3/2}i^{3/2}}{2n}.\end{eqnarray}

Similarly, when $c\leq 1,$ without making an assumption on whether $ci\leq 1$ or not,

(5.4) \begin{eqnarray}\mathbb{P}\bigl(\widetilde X_{k+1}\neq \widetilde Z_{k+1}\bigl|\widetilde X_k=\widetilde Z_k=i\bigr)\leq \frac{ci}{2n}+\frac{ci^2}{2n}+\frac{c^2i^2}{2n}\leq \frac{3ci^2}{2n}.\end{eqnarray}

Recall $\tau=\inf\bigl\{k\in\mathbb{N}: \widetilde X_k\neq \widetilde Z_k\bigr\}$ and let

\begin{eqnarray*}\sigma=\inf\bigl\{k\in \mathbb{N}: \widetilde Z_k=0\bigr\}.\end{eqnarray*}

To evaluate the distribution function of T we will use the following inequalities:

(5.5) \begin{align}\nonumber\mathbb{P}(T\leq m)&\leq \mathbb{P}(T\leq m,\tau>m)+\mathbb{P}(\tau\leq m)\\[3pt]&\leq\mathbb{P}(\sigma\leq m)+\mathbb{P}(\tau\leq m)\end{align}

and

(5.6) \begin{eqnarray}\mathbb{P}(T\leq m)\geq (\sigma\leq m).\end{eqnarray}

The latter inequality holds because the avalanche process X is stochastically dominated by the branching process $Z^{(c)}$, by virtue of Proposition 3.1.

The following lemma summarizes results of [Reference Agresti1] regarding the distribution (subdistribution if the process is supercritical) function $\mathbb{P}(\sigma \leq m).$ Specific bounds for the extinction time of a Poisson distribution are derived in Theorem 2 of [Reference Agresti1].

Recall $\alpha_c$ from (4.4) and s(c),r(c) from (5.1).

We next estimate $\mathbb{P}(\tau\leq m).$ For $k\in\mathbb{N}$ let $W_k=\max_{0\leq i \leq k} \widetilde Z_i.$ By the Markov property of $\bigl(\widetilde X_k,\widetilde Z_k\bigr)_{k\in \mathbb{Z}_+},$ for any $J\in\mathbb{N}$ we have

(5.7) \begin{align}\nonumber\mathbb{P}(\tau\leq m)&\leq \mathbb{P}(\tau\leq m\ \textrm{and}\ W_m< J)+\mathbb{P}(W_m\geq J)\\\nonumber&\leq\mathbb{P}\Bigl(\bigcup_{k=1}^m \bigl\{\widetilde X_{k-1}=\widetilde Z_{k-1}<J,\,\widetilde X_k\neq \widetilde Z_k\bigr\}\Bigr)+\mathbb{P}(W_m\geq J)\\&\leq m\cdot \mathbb{P}\bigl(\widetilde X_k\neq \widetilde Z_k\,\bigr|\,\widetilde X_{k-1}=\widetilde Z_{k-1}<J\bigr) +\mathbb{P}(W_m\geq J).\end{align}

The first part of the next lemma is an improved version of Lemma 5.6 in [Reference Chumley, Aydogmus, Matzavinos and Roitershtein13].

Proof of Lemma 5.3.

1. (i) Let $U_k\,{:\!=}\,c^{-k}\widetilde Z_k.$ Then $(U_k)_{k\in\mathbb{Z}_+}$ is a martingale with respect to its natural filtration. For any $\theta>0,$ $f(x)=e^{\theta x}$ is a convex function, and hence the sequence $e^{\theta U_k},$ $k\in \mathbb{Z}_+,$ forms a submartingale. Hence, by Doob’s maximal inequality,

   \begin{align*}\mathbb{P}(W_m\geq xc^m)\leq \bigl(\max_{0\leq k \leq m} e^{\theta U_k}\geq e^{\theta x}\bigr)\leq e^{-\theta x}\mathbb{E}\bigl(e^{\theta U_m}\bigr).\end{align*}
   The result now follows from Theorem 4 in [Reference Athreya4], which states that $\sup_{m\in\mathbb{N}} \mathbb{E}\bigl(e^{\theta U_m}\bigr)<+\infty$ for some $\theta>0.$

   (ii) and (iii) These are direct implications of Theorem 2 in [Reference Lindvall27].

The claims in the theorem now follow from combining the bounds in (5.3)–(5.7) with the estimates in Lemmas 5.2 and 5.3. In the case when $c=1,$ the optimization problem over the optimal choice of the parameter J in the upper bound for $\mathbb{P}(T\leq m)$ can be solved explicitly, and we choose

\begin{align*}J=\left(\frac{n i_0^2}{3m}\right)^{1/3}.\\[-35pt]\end{align*}

We will next consider a family of avalanches under Assumption 3.1. For $J>0$ let

(5.8) \begin{eqnarray}T^{(n)}=\min\bigl\{k\in\mathbb{N}:X^{(n)}_k=0\bigr\}.\end{eqnarray}

The following result shows that the bounds in Lemma 5.2 hold asymptotically for the avalanche model.

Proof. Informally speaking, the obvious strategy for the proof is to take n to infinity in the conclusions of Theorem 5.1. Let $\alpha_1=r(1)=s(1)=1.$ Observe that for any $\lambda>0,$ $\lim_{\mu\to \lambda} \alpha_\mu=\alpha_\lambda,$ and for any $\lambda \in [0,1],$ $\lim_{\mu\to \lambda} s(\mu)=s(\lambda)$ and $\lim_{\mu\to \lambda} r(\mu)=r(\lambda),$ where in the case $\lambda=1$ the latter limits are understood as left limits. Furthermore, the limits in Lemma 5.2 are continuous functions of c. This is evident for $c\neq 1,$ and for $c=1$ by L’Hospital’s rule we have

\begin{eqnarray*}\lim_{c\to 1^-} \frac{s(c)(1-c^m)}{s(c)-c^m}=\lim_{c\to 1^-}\frac{-(1-c^m)-m(2-c)c^{m-1}}{-1-(m+1)c^m}=\frac{m}{m+2}\end{eqnarray*}

and

\begin{eqnarray*}\lim_{c\to 1^-} \frac{r(c)(1-c^m)}{r(c)-c^m}&=&\lim_{c\to 1^-}\frac{(e^{-c}-ce^{-c})(1-c^m)-me^{-c}c^m}{e^{-c}-ce^{-c}-mc^{m-1}(e^{-c}-1+c)-c^m(-e^{-c}+1)}\\&=&\frac{me^{-1}}{me^{-1}-e^{-1}+1}=\frac{m}{m+e-1},\end{eqnarray*}

as desired.

For $k\in\mathbb{N}$ and $\mu>0,$ let $W^{(\mu)}_k=\max_{0\leq i \leq k} Z^{(\mu)}_i.$ In comparison with the proof of Theorem 5.1, the only nuance here is that before letting n go to infinity, the term $\mathbb{P}\bigl(W^{(np_n)}_m\geq J\bigr)$ in (5.7), which a priori is not uniform in n, should be replaced with $\mathbb{P}\bigl(W^{(c)}_m\geq J_n\bigr),$ where c is the constant introduced in Condition 3.1 and $J_n$ is a suitable sequence which in particular satisfies $\lim_{n\to\infty}J_n=+\infty.$ We remark that technically, this part of the proof makes the claim a corollary to a slight modification of Theorem 5.1 rather than to the theorem itself.

Another interesting consequence of Theorem 5.1 is the following result, which shows that on the right scale, the asymptotic behavior of the avalanche model coincides with that of the branching process (see, for instance, the compact introduction section in [Reference Lyons, Pemantle and Peres29] for the underlying branching process estimates). More precisely, the next corollary evaluates the asymptotic behavior of $\mathbb{P}\bigl(T^{(n)}> m_n\bigr)$ for a range of ‘moderately fast-growing’ sequences $m_n.$

Proof. In each of the three cases under consideration (namely $\lambda>1$, $\lambda<1$, $\lambda=1$), the sequence $\beta_n$ is chosen in such a way that the term corresponding to $\mathbb{P}(W_n\geq m)$ in the key estimate (5.7) and consequently in the conclusions of Theorem 5.1 is dominated by lower and upper bounds for $\mathbb{P}(\sigma\leq m)$ which are supplied by Lemma 5.2. As a result, the term vanishes as n goes to infinity, and the asymptotic behavior of $\mathbb{P}\bigl(T^{(n)}>m_n\bigr)$ coincides with that of its counterpart for branching processes. The only subtlety is in the proof of Part (iii), where the asymptotic behavior exhibited by the lower and upper bounds in Lemma 5.2 do not match:

\begin{equation*}\lim_{n\to\infty} m\left\{1-\left(\frac{m}{m+2}\right)^{i_0}\right\}=2i_0\end{equation*}

while

\begin{eqnarray*}\lim_{n\to\infty} m\left\{1-\left(\frac{m}{m+e-1}\right)^{i_0}\right\}=(e-1)i_0.\end{eqnarray*}

However, Kolmogorov’s estimate (see, for instance, Theorem C in [Reference Lyons, Pemantle and Peres29]) ensures that $\lim_{m\to\infty} \mathbb{P}(\sigma>m)=2i_0.$ In view of (5.7) this yields the desired result in the case $\lambda=1.$ The proof of the corollary is complete.

6. Deterministic approximation

The branching approximation reflects the dynamics of the process suitably when the latter is subcritical and the size of avalanche decays exponentially fast, similarly to the branching approximation process. In contrast, in the supercritical regime, namely when Condition 2.1 is satisfied with $c>1$ or Assumption 3.1 holds with $\lambda>1,$ the results in Section 2 tend to provide little or no information on the asymptotic behavior of the model, and even a partial improvement of this situation is clearly desirable. However, even in the supercritical case, one should expect that with $\varepsilon \ll 1,$ the branching approximation is adequate for $\sim \varepsilon \log n$ of first steps while the ratio $\frac{X_k}{n}$ remains small; cf. Section 6.3.1 of [Reference Durrett16]. On the other hand, once the ratio becomes of order one, one can expect that the dynamics of the model will be almost deterministic and will follow a ‘mean-field’ equation due to the law of large numbers. This section is devoted to the study of this deterministic approximation of the model. The relevance of this regime to a supercritical avalanche model is formally elucidated by the results in Section 4.

The rest of the section is divided into three subsections. Subsection 6.1 discusses a heuristic computation based on a complete decoupling of a weak dependence between the network nodes, which sheds additional light on the nature of the deterministic approximation. In Subsection 6.2 the deterministic dynamical system is formally introduced and studied. Interestingly enough, the underlying dynamics resembles but is richer than that of the logistic equation. This was already observed in [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15, Reference Longini28]. Subsection 6.3 addresses the question of how long the trajectories of an avalanche Markov chain and the corresponding deterministic dynamical system remain close each to other, provided they began at the same point.

6.1. Extra motivation and a precise comparison result

Our motivation in this section comes partially from the following heuristic computation, which is an adaptation to our framework of the one proposed in [Reference Gleeson18] for a model of avalanches in a Boolean network. Consider an avalanche Markov chain $(X_k)_{k\in\mathbb{Z}_+}$ with transition kernel introduced in (1.2), and recall $\mathcal{E}_k$ from (2.3) and $\xi_k(x)$ from (2.6). If the random variables $\{\mathcal{E}_k(x):x\in V_n\}$ were i.i.d. for some $k\in\mathbb{Z}_+,$ we could relate their common value $\xi_k$ to the common value $\xi_{k+1}$ of $\xi_{k+1}(x),$ $x\in V_n,$ by means of the following recursion [Reference Gleeson18]:

(6.1) \begin{eqnarray}\nonumber\xi_{k+1}&=&\bigl(1-\xi_k\bigr)\sum_{i=0}^{n-1}\binom{n-1}{i}\xi_k^i(1-\xi_k)^{n-1-i}\bigl(1-q^{i}\bigr)\\&=&(1-\xi_k)\bigl[1-(q\xi_k+1-\xi_k)^n\bigr]=(1-\xi_k)\bigl[1-(1-p\xi_k)^n\bigr].\end{eqnarray}

Consequently, the following inequality would hold:

(6.2) \begin{eqnarray}\xi_{k+1}\geq (1-\xi_k)\bigl(1-e^{-pn\xi_k}\bigr)= (1-\xi_k)\bigl(1-e^{-c\xi_k}\bigr),\end{eqnarray}

where we write $c=np.$ We remark that similar decoupling arguments are often used in the physics literature to justify a mean-field approximation in a complex locally tree-like network (hence a fairly weak dependence between the nodes); see, for instance, [Reference Larremore, Carpenter, Ott and Restrepo24–Reference Larremore, Shew, Ott and Restrepo26, Reference Rämö, Kesseli and Yli-Harja33, Reference Rämö, Kauffman, Kesseli and Yli-Harja34].

Of course, the $\mathcal{E}_k(x)$ are not independent and in general are not identically distributed (the latter depends on whether the distribution of $X_0$ is exchangeable or not; cf. Proposition 2.3). For a different but somewhat related model of avalanches, it is argued in [Reference Gleeson18] (see the comment [Reference Athreya4] in the references section there) that the above i.i.d. assumption ‘is true for large networks with well-behaved degree distributions, but excludes networks with hubs which output to a substantial fraction of the nodes.’ It is of interest to note that, in agreement with this general heuristic principle, a suitable modification of (6.1) and, consequently, of (6.2) serves in a certain rigorous sense as a good approximation for the avalanche model. This is accomplished in (6.11) below.

We will now proceed with a rigorous modification of the above heuristic calculation. For $\alpha>0$ let $g_\alpha(x)=(1-x)(1-e^{-\alpha x}),$ $x\in [0,1],$ and

\begin{eqnarray*}\chi_\alpha=\max_{x\in(0,1)}g_\alpha(x)\quad \textrm{and}\quad\nu_\alpha=\underset{x\in(0,1)}{\textrm{argmax}}\ g_\alpha(x).\end{eqnarray*}

Note that $\nu_\alpha$ is uniquely defined for all $\alpha>0.$ Further, if $\alpha>1,$ let $\zeta_\alpha$ be the unique solution on (0,1) to the fixed point equation $g_\alpha(\zeta_\alpha)=\zeta_\alpha.$

Notice that $\zeta_\alpha=g_\alpha(\zeta_\alpha)<1-\zeta_\alpha$ implies that $\zeta_\alpha<1/2.$ The following lemma summarizes other basic properties of the function $g_\alpha$ that we are going to use (see also Theorem 6.3 in Subsection 6.2 below). The proof of the lemma is omitted; Fig. 3 below illustrates the dependence of the graph of $g_\alpha$ on the parameter $\alpha$, as well as basic properties of $g_\alpha$. For more details see, for instance, [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15], where the function $g_\alpha$ is studied systematically in a similar context.

Figure 3.

Graph of the function $g_\alpha(x)=(1-x)(1-e^{-\alpha x})$ for several values of the parameter $\alpha,$ including the critical branching value $\alpha=1,$ for which $g_1'(0)=1,$ and $\alpha=2.46742$, which is a close approximation to the transitional value $\alpha_{\rm tr}$.

Lemma 6.1.

1. (i) The function $g_\alpha(x)$ is increasing on $(0,\nu_\alpha)$ and decreasing on $(\nu_\alpha,1).$

2. (ii) If $\alpha\leq 1,$ then $g_\alpha(x)<x$ for all $x\in (0,1).$

3. (iii) If $\alpha>1,$ then $g_\alpha(x)>x$ on $(0,\zeta_\alpha)$ and $g_\alpha(x)<x$ on $(\zeta_\alpha,1).$

4. (iv) We have $g_\alpha(x)<\alpha x$ for all $\alpha>0$ and $x\in (0,1).$

5. (v) There exists a transitional value $\alpha_{\rm tr}>1$ such that $\nu_\alpha>\zeta_\alpha$ if and only if $\alpha<\alpha_{\rm tr}.$

Numerical simulations indicate that $\alpha_{\rm tr}\approx 2.46742.$ The role of the threshold in the dynamics of the model is discussed in more detail in [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15, Reference Rastegar and Roitershtein35].

For an avalanche model with transition kernel (1.2) and fixed network size $n\in\mathbb{N},$ let

(6.3) \begin{eqnarray}\varphi_k=\frac{1}{n}\mathbb{E}(X_k), \qquad k\in\mathbb{Z}_+.\end{eqnarray}

Note that if the distribution of $X_0$ is exchangeable (i.e. invariant with respect to permutation of the network nodes), then $\varphi_k=\xi_k(x),$ where $\xi_k$ is defined in (2.6), for all $k\in\mathbb{Z}_+$ and $x\in V_n.$ By the bounded convergence theorem, for any fixed $n\in\mathbb{N},$ $\lim_{k\to\infty} \varphi_k=0$ regardless of the choice of parameter $p>0.$ We have the following theorem.

Theorem 6.1. Let $(X_k)_{k\in\mathbb{Z}_+}$ be an avalanche model with transition kernel (1.2). Define

(6.4) \begin{eqnarray}\alpha=-n\log(1-p).\end{eqnarray}

The following hold:

1. (i) Let $(\phi_k)_{k\in\mathbb{Z}_+}$ be defined recursively by setting $\phi_0=\frac{1}{n}\mathbb{E}(X_0)$ and

   \begin{eqnarray*}\phi_{k+1}=1-e^{-\alpha \phi_k}.\end{eqnarray*}
   Then $\varphi_k\leq \phi_k$ for all $k\in\mathbb{Z}_+.$

2. (ii) We have $\sup_{k\in\mathbb{Z}_+} \varphi_k\leq \max\bigl\{\phi_0,\chi_\alpha\bigr\}.$

3. (iii) Let $(\psi_k)_{k\in\mathbb{Z}_+}$ be defined recursively by setting $\psi_0=\frac{1}{n}\mathbb{E}(X_0)$ and

   \begin{eqnarray*}\psi_{k+1}=g_\alpha (\psi_k).\end{eqnarray*}
   If $\alpha\leq \alpha_{\rm tr}$ and $\varphi_0\leq \zeta_\alpha,$ then
   \begin{eqnarray*}\varphi_k\leq \psi_k \quad \forall\,k\in\mathbb{Z}_+, \qquad \textrm{and moreover,}\qquad\sup_{k\in\mathbb{Z}_+} \varphi_k\leq \zeta_\alpha \leq \chi_\alpha.\end{eqnarray*}

Note that in view of the fact that $1-e^{-\alpha \phi_k}<\alpha \phi_k$ and the result in Part (iv) of Lemma 6.1, the theorem improves the result in the basic Corollary 2.1 even in the subcritical case $\alpha<1.$ The relevance of the phase transition at $\alpha_{\rm tr}$ to the dynamics of the avalanche Markov chain is further discussed in Section 6.2 below, in the paragraph following Corollary 6.1.

Proof of Theorem 6.1.

1. (i) It follows from (1.2), (6.4), and Jensen’s inequality (applied in the last step to the function $g(x)=(1-x)(1-e^{-\alpha x})$, which is concave on (0,1)) that

   (6.5) \begin{eqnarray}\nonumber\varphi_{k+1}(x)&=&\mathbb{E}\Bigl[\Bigl(1-\frac{X_k}{n}\Bigr)\bigl(1-q^{X_k}\bigr)\Bigr]=\mathbb{E}\Bigl[\Bigl(1-\frac{X_k}{n}\Bigr)\Bigl(1-\bigl(1-p\bigr)^{X_k}\Bigr)\Bigr]\\[3pt] \nonumber&=&\mathbb{E}\left[\left(1-\frac{X_k}{n}\right)\left(1-\Bigl\{\bigl(1-p\bigr)^{p}\Bigr\}^{\frac{X_k}{p}}\right)\right]\\[3pt] &=&\mathbb{E}\Bigl[\Bigl(1-\frac{X_k}{n}\Bigr)\Bigl(1-e^{-\frac{\alpha X_k}{n}}\Bigr)\Bigr]\leq (1-\varphi_k)(1-e^{-\alpha \varphi_k}).\end{eqnarray}
   In particular, $\varphi_{k+1}\leq 1-e^{-\alpha \varphi_k}.$ Since $f(x)=1-e^{-\alpha x}$ is a monotone increasing function, an induction argument shows that $\varphi_k\leq \phi_k$ for all $k\in\mathbb{Z}_+.$

2. (ii) The claim is immediate from (6.5).

3. (iii) Observe that $g_\alpha$ is monotone increasing on $(0,\zeta_\alpha)$ for $\alpha\leq \alpha_{\rm tr}.$ Therefore it follows from (6.5), via induction on k, that $\varphi_k\leq \psi_k\leq \zeta_\alpha$ for all $k\in\mathbb{Z}_+.$ Moreover, if $\alpha\leq 1$ then $\psi(x)<x$ on $(0,\zeta_\alpha),$ and hence $\psi_{k+1}\leq \psi_k$ for all $k\in\mathbb{Z}_+,$ while if $\alpha\in (1,\alpha_{\rm tr}]$ then $\psi(x)>x$ on $(0,\zeta_\alpha),$ and hence $\psi_{k+1}\geq \psi_k$ for all $k\in\mathbb{Z}_+.$ Since $g_\alpha$ has either one fixed point, at zero (for $\alpha\leq 1$), or two, at zero and $\zeta_\alpha$ (for $\alpha>1$), this implies that $\lim_{k\to\infty} \varphi_k=0$ for $\alpha\leq 1$, while

   \begin{equation*}\sup_{k\in\mathbb{Z}_+} \varphi_k\leq \sup_{k\in\mathbb{Z}_+} \psi_k\leq \lim_{k\to\infty} \varphi_k=\zeta_\alpha\end{equation*}
   in the case that $\alpha \in (1,\alpha_{\rm tr}].$

6.2. Deterministic approximation for large networks

We turn now to a study of an ensemble of avalanche models which satisfies Assumption 3.1. First we discuss a direct implication of the results in Theorem 6.1 for the asymptotic behavior of the expected fraction of excited nodes in a network of the ensemble. The main results of this section are stated afterwards in Theorems 6.2 (regarding asymptotic behavior of $\frac{1}{n}X^{(n)}_k$ for large n) and Corollary 6.2 (a consequence of the results in Theorem 6.2 for the heterogeneity of a network in the ensemble).

First, observe that the continuity of the results in Theorem 6.1 in the parameters $\alpha$ and the initial data $\varphi_0$, together with the monotonicity of $g_\alpha$ on the interval $\bigl(0,\zeta_\alpha\bigr)$, leads to the following corollary to the theorem.

Corollary 6.1. Let Assumption 3.1 hold, and suppose that the following limit exists and belongs to $(0,1):$

\begin{eqnarray*}\varepsilon_0=\lim_{n\to\infty} \frac{1}{n}\mathbb{E}\bigl(X^{(n)}_0\bigr).\end{eqnarray*}

For $n\in\mathbb{N}$ and $k\in\mathbb{Z}_+$ let

(6.6) \begin{eqnarray}\varphi_{n,k}=\frac{1}{n}\mathbb{E}\bigl(X^{(n)}_k\bigr)\end{eqnarray}

be the counterpart of $\varphi_k$ introduced for a single network in (6.3).

Then the following hold:

1. (i) Let $(\phi_k)_{k\in\mathbb{Z}_+}$ be defined recursively by setting $\phi_0=\varepsilon_0$ and

   \begin{eqnarray*}\phi_{k+1}=1-e^{-\lambda \phi_k}.\end{eqnarray*}
   Then $\limsup_{n\to\infty}\varphi_{n,k}\leq \phi_k$ for all $k\in\mathbb{Z}_+.$

2. (ii) We have $\limsup_{n\to\infty} \varphi_{n,k}\leq \chi_\lambda$ for all $k\in\mathbb{N}.$

3. (iii) Recall $\alpha_{\rm tr}$ from Lemma 6.1. Let $(\psi_k)_{k\in\mathbb{Z}_+}$ be defined recursively by setting

   (6.7) \begin{eqnarray}\psi_0=\varepsilon_0\qquad \textrm{and}\qquad \psi_{k+1}=g_\lambda (\psi_k).\end{eqnarray}
   If $\lambda<\alpha_{\rm tr}$ and $\varepsilon_0<\zeta_\lambda,$ then
   \begin{eqnarray*}\limsup_{n\to\infty}\varphi_{n,k}\leq \psi_k\qquad \textrm{and}\qquad\limsup_{n\to\infty} \varphi_{n,k}\leq \zeta_\lambda \leq \chi_\lambda\end{eqnarray*}
   for all $k\in\mathbb{Z}_+.$

Let

(6.8) \begin{eqnarray}x_{n,k}=\frac{X^{(n)}_k}{n},\qquad n\in\mathbb{N},\,k\in\mathbb{Z}_+.\end{eqnarray}

Under Assumption 3.1, define the asymptotic branching factor as a function $b:(0,1)\to\mathbb{R}$ by setting

\begin{eqnarray*}b(x)=\lim_{n\to\infty} \frac{1}{x}\mathbb{E}(x_{n,k+1}\,|\,x_{n,k}=x)=\frac{g_\lambda(x)}{x},\qquad x\in (0,1).\end{eqnarray*}

It turns out (see [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15] for details) that the behavior of the sequence $b(\psi_k)$ for $\alpha\in(1,\alpha_{\rm tr})$ differs qualitatively from the behavior for $\alpha>\alpha_{\rm tr}$; that is, a secondary phase transition in the avalanche model occurs at the transitional value $\alpha_{\rm tr}.$ For instance, if $\alpha\in (1,\alpha_{\rm tr})$ and $\psi_0$ is sufficiently small, then $b(\psi_k)$ increases monotonically to one as $k\to\infty.$ In contrast, this is not necessarily true when $\alpha>\alpha_{\rm tr}.$ For example, in the case $\alpha>\alpha_{\rm tr},$ if for some $m\in\mathbb{N},$ $\psi_0<\nu_\alpha$ and $\psi_m=\nu_m,$ then the sequence $\psi_k$ increases monotonically at the first m steps and then converges to one by oscillating between values which are larger and smaller than one [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15, p. 81]. Note that by choosing m sufficiently large, we can place $\psi_0$ as close to zero (the only fixed point of the equation $x=g_\alpha(x)$ other than $\zeta_\alpha$) as we wish.

Let $\overset{P}{\to}$ denote convergence in probability as the size of the network n goes to infinity. The following theorem is an adaptation to our setup of Theorems 1 and 3 in [Reference Buckley and Pollett11] (see also [Reference Klebaner and Nerman23] for earlier similar results).

Theorem 6.2. Let Assumption 3.1 hold. Recall $x_{n,k}$ from (6.8) and suppose that

(6.9) \begin{eqnarray}x_{n,0}\overset{P}{\to}\psi_0\end{eqnarray}

for some constant $\psi_0\in (0,1).$ Then the following hold:

1. (a) We have $x_{n,k}\overset{P}{\to} \psi_k$ for all $k\in\mathbb{Z}_+,$ where $\psi_{k+1}=g_\lambda(\psi_k).$

2. (b) Let

   \begin{eqnarray*}y_{n,k}=\sqrt{n}(x_{n,k}-\psi_k),\qquad n\in\mathbb{N},\,k\in\mathbb{Z}_+,\end{eqnarray*}
   and set
   \begin{eqnarray*}v(x)\,{:\!=}\,g_\lambda(x)e^{-\lambda x}=(1-x)e^{-\lambda x}(1-e^{-\lambda x}),\qquad x\in [0,1].\end{eqnarray*}

Suppose that in addition to (6.9), $y_{n,0}$ converges weakly, as n goes to infinity, to some (possibly random) $Y_0.$ Then the sequence $y^{(n)}\,{:\!=}\,(y_{n,k})_{k\in\mathbb{Z}_+}$ converges in distribution, as n goes to infinity, to a time-inhomogeneous Gaussian AR(1) sequence $(Y_k)_{k\in\mathbb{Z}_+}$ defined by

(6.10) \begin{eqnarray}Y_{k+1}=g^{\prime}_\lambda(\psi_k)Y_k+e_k,\end{eqnarray}

where the $e_k,$ $k\in\mathbb{Z}_+,$ are independent Gaussian variables, with each $e_k$ distributed as $N\bigl(0,v(\psi_k)\bigr).$

Proof. Let ${\mathcal Z}\,{:\!=}\,\bigl\{{\mathcal Z}_{k,i}^{(n,j)}:n,i,j\in\mathbb{N},k\in\mathbb{Z}_{+}\bigr\}$ be a collection of independent Bernoulli variables with

\begin{eqnarray*}\mathbb{P}\bigl({\mathcal Z}_{k,i}^{(n,j)}=1\bigr)=1-q_n^j\qquad \textrm{and}\qquad \mathbb{P}\bigl({\mathcal Z}_{k,i}^{(n,j)}=0\bigr)=q_n^j.\end{eqnarray*}

Thus, without loss of generality, we can assume that

\begin{eqnarray*}X^{(n)}_{k+1}=\sum_{i=1}^{n-X^{(n)}_k} {\mathcal Z}_{k,i}^{(n,X^{(n)}_k)}.\end{eqnarray*}

Let ${\mathcal Y}\,{:\!=}\,\bigl\{{\mathcal Y}_{k,i}^{(n,j)}:n,i,j\in\mathbb{N},k\in\mathbb{Z}_+\bigr\}$ be another collection of independent Bernoulli variables defined on the same probability space, and such that

\begin{eqnarray*}\mathbb{P}\bigl({\mathcal Y}_{k,i}^{(n,j)}=1\bigr)=1-e^{-\frac{\lambda j}{n}}\qquad \textrm{and}\qquad \mathbb{P}\bigl({\mathcal Y}_{k,i}^{(n,j)}=0\bigr)=e^{-\frac{\lambda j}{n}}.\end{eqnarray*}

For $n\in\mathbb{N},$ let $\alpha_n=np_n.$ By using the maximal coupling for two Bernoulli variables, we can and will assume that the pairs $\bigl({\mathcal Z}_{k,i}^{(n,j)},{\mathcal Y}_{k,i}^{(n,j)}\bigr)$ are independent $\{0,1\}^2$-random variables, and

\begin{eqnarray*}\mathbb{P}\bigl({\mathcal Z}_{k,i}^{(n,j)}\neq {\mathcal Y}_{k,i}^{(n,j)}\bigr)=\Bigl|e^{-\frac{\lambda j}{ n}}-\Bigl(1-\frac{\alpha_n}{n}\Bigr)^j\Bigr|.\end{eqnarray*}

For $n\in\mathbb{N},$ define a new sequence ${\widetilde X}^{(n)}=\bigl({\widetilde X}^{(n)}_k\bigr)_{k\in\mathbb{Z}_+}$ by setting ${\widetilde X}^{(n)}_0=\psi_0$ and

\begin{eqnarray*}{\widetilde X}^{(n)}_{k+1}=\sum_{i=1}^{n-{\widetilde X}^{(n)}_k} {\mathcal Y}_{k,i}^{(n,{\widetilde X}^{(n)}_k)}.\end{eqnarray*}

Theorems 1 and 3 in [Reference Buckley and Pollett11] ensure that the results in the theorem, both the law of large numbers and the central limit theorem, hold if we replace $X^{(n)}$ with ${\widetilde X}^{(n)}.$ Thus, in order to prove the theorem, it suffices to show that $\bigl(X^{(n)}_k-{\widetilde X}_k^{(n)}\bigr)\overset{P}{\to} 0$ for all $k\in \mathbb{Z}_+$ (see, for instance, Remark (i) on p. 60 of [Reference Buckley and Pollett11] and/or the last paragraph in the proof of Theorem 3 there, which both assert that the main result of [Reference Karr20] goes through to a non-homogeneous chain setting, and therefore the weak convergence in question is implied by the convergence of transition kernels). To this end, observe that

\begin{align*}\mathbb{E}\left( \frac{1}{n}\Bigl|X^{(n)}_{k+1}-{\widetilde X}_{k+1}^{(n)}\Bigr|\right)&\leq\frac{1}{n}\mathbb{E}\left(\sum_{i=1}^{n-X^{(n)}_k} \Bigl|{\mathcal Z}_{k,i}^{(n,X^{(n)}_k)}-{\mathcal Y}_{k,i}^{(n,X^{(n)}_k)}\Bigr|\right)\\ &\leq \mathbb{E}\Bigl(\Bigl|{\mathcal Z}_{k,1}^{(n,X^{(n)}_k)}-{\mathcal Y}_{k,1}^{(n,X^{(n)}_k)}\Bigr|\Bigr)=\mathbb{P}\Bigl({\mathcal Z}_{k,1}^{(n,X^{(n)}_k)}\neq {\mathcal Y}_{k,1}^{(n,X^{(n)}_k)}\Bigr)\\ &= \mathbb{E}\Bigl(\Bigl|e^{-\frac{\lambda X^{(n)}_k}{ n}}-\Bigl(1-\frac{\alpha_n}{n}\Bigr)^{X^{(n)}_k}\Bigr|\Bigr),\end{align*}

and hence the claim can be proved by induction, using the bounded convergence theorem.

Part (a) of the above result and the bounded convergence theorem imply that

(6.11) \begin{eqnarray}\lim_{n\to\infty} \varphi_{n,k}=\psi_k\qquad \forall\,k\in\mathbb{Z}_+,\end{eqnarray}

where $\varphi_{n,k}$ is defined in (6.3). Note that this limit identity is reminiscent of the heuristic (6.1) and (6.2) in our framework.

It is not hard to prove (cf. Remark (v) on p. 61 of [Reference Buckley and Pollett11]; see also [Reference Klebaner and Nerman23]) that if $\psi_0\overset{P}{\to} \zeta_\lambda$ in the statement of Theorem 6.2 and, in addition, $y_{n,0}$ converges weakly, as n goes to infinity, to some $Y_0,$ then the linear recursion (6.10) can be replaced with

\begin{eqnarray*}Y_{k+1}=g^{\prime}_\lambda(\zeta_\lambda)Y_k+\widetilde e_k,\end{eqnarray*}

where $\widetilde e_k,$ $k\in\mathbb{Z}_+,$ are i.i.d. Gaussian variables, with each $\widetilde e_k$ distributed as $N\bigl(0,v(\zeta_\lambda)\bigr).$ One then can show (see the proof of Theorem 6.4 below) that $|g_\lambda^{\prime}(\zeta_\lambda)|<1,$ and hence the Markov chain $Y_k$ has a stationary distribution; see [Reference Buckley and Pollett11, p. 61] for more details.

Recall $H_k$ from (2.2) and define a normalized heterogeneity $h_{n,k}$ by

\begin{eqnarray*}h_{n,k}=\frac{H_{n,k}}{2n(n-1)},\qquad n\in\mathbb{N},k\in\mathbb{Z}_+.\end{eqnarray*}

Notice that if the distribution of $X^{(n)}_0$ is exchangeable, then $x_{n,k}$ is a probability that two nodes in the network generating $X^{(n)}_k$ randomly chosen at time k have different types. The following is immediate from Theorem 6.2.

Corollary 6.2. Under the conditions of Theorem 6.2, the following hold:

1. (a) For $x\in [0,1]$ let $r(x)=\frac{1}{2}x(1-x).$ Then $h_{n,k}\overset{P}{\to} r(\psi_k)$ for all $k\in\mathbb{Z}_+.$

2. (b) Let

   \begin{eqnarray*}{\widetilde y}_{n,k}=\sqrt{n}\bigl(h_{n,k}-r(\psi_k)\bigr),\qquad n\in\mathbb{N},\,k\in\mathbb{Z}_+.\end{eqnarray*}

Suppose that in addition to (6.8), $y_{n,0}=\sqrt{n}\bigl(x_{n,0}-\psi_0\bigr)$ converges weakly, as n goes to infinity, to some (possibly random) $Y_0.$ Then the sequence ${\widetilde y}^{(n)}\,{:\!=}\,({\widetilde y}_{n,k})_{k\in\mathbb{Z}_+}$ converges in distribution, as n goes to infinity, to a time-inhomogeneous Gaussian AR(1) sequence $(\widetilde Y_k)_{k\in\mathbb{Z}_+},$ where $\widetilde Y_k=\frac{1}{2}(1-2\psi_k)Y_k$ and $Y_k$ is defined in (6.10).

Proof. The convergence in probability of $h_{n,k}$ follows from the continuous mapping theorem and the result in Part (a) of Theorem 6.2. To prove the weak convergence of ${\widetilde y}_{n,k},$ write

\begin{eqnarray*}{\widetilde y}_{n,k}=\sqrt{n}\bigl(r(x_{n,k})-r(\psi_k)\bigr)+\frac{\sqrt{n}}{2}X_{n,k}(n-X_{n,k})\Bigl(\frac{1}{n(n-1)}-\frac{1}{n^2}\Bigr)\end{eqnarray*}

and observe that

\begin{eqnarray*}\Bigl|\frac{\sqrt{n}}{2}X_{n,k}(n-X_{n,k})\Bigl(\frac{1}{n(n-1)}-\frac{1}{n^2}\Bigr)\Bigr|\leq \frac{\sqrt{n}}{2(n-1)} \longrightarrow_{n\to\infty} 0.\end{eqnarray*}

Furthermore, by the mean value theorem,

\begin{eqnarray*}\sqrt{n}\bigl(r(x_{n,k})-r(\psi_k)\bigr)=r'(x^*_{n,k})\sqrt{n}(x_{n,k}-\psi_k)=(1-2x_{n,k})\sqrt{n}(x_{n,k}-\psi_k)\end{eqnarray*}

for some $x^*_{n,k}$ between $x_{n,k}$ and $\psi_k.$ In view of the result in Part (b) of Theorem 6.2, the claim in Part (b) of this theorem now follows from another application of the continuous mapping theorem.

The discrete-time dynamical system $\psi_k$ is studied in detail in [Reference Cooke, Calef, Level, Lakshmikantham, Press and York15]. In particular, the following result is proved there.

Theorem 6.3. ([Reference Gleeson and Durrett15, Theorem 1]) Let $\lambda>0$ and $\varepsilon_0\in(0,1)$ be given, and let the sequence $(\psi_k)_{k\in\mathbb{Z}_+}$ be defined as in (6.7). Then the following hold

1. (i) If $\lambda\in (0,1],$ then $\lim_{k\to\infty} \psi_k=0$ for any $\psi_0\in [0,1].$

2. (ii) If $\lambda>1,$ then for all $\psi_0\in [0,1],$ $\lim_{k\to\infty} \psi_k=\zeta_\lambda,$ where $\zeta_\lambda\in (0,1/2)$ is the unique positive solution to the fixed point equation $g_\lambda(x)=x.$

Together with the results in Section 4 and Theorem 6.2, this theorem implies that when n is large, with high probability, a supercritical Markov chain $X^{(n)}$ will eventually be trapped for a long time in a neighborhood of the global stable point $\zeta_\lambda$ of the map $g_\lambda.$ The next section is devoted to the proof of a certain qualitative form of this informal observation.

6.3. Comparison of the stochastic and deterministic trajectories

The main results of this section are stated in Theorem 6.4 (supercritical case) and Theorem 6.5 (critical and subcritical cases).

Theorem 6.2 suggests that when both n and the size of the first generation $X^{(n)}_0$ in the avalanche process are substantially large, the trajectory of the deterministic sequence $\psi_k$ can serve as a good approximation to the path of the Markov chain $X^{(n)}_k.$ Note however that, at least for a supercritical process, two trajectories cannot in principle stay close each to other forever, since while the latter converges to zero with probability one, the former tends to a non-zero limit by virtue of Theorem 6.3.

The following theorem offers some insight into the duration of the time when the deterministic and the stochastic paths stay fairly close to each other, before they first become significantly separated from each other. The theorem is a suitable modification of some results in [Reference Aydogmus6]. In words, the theorem asserts that if Assumption 3.1 holds with $\lambda>1$, if the scaled initial population $x_{n,0}$ is close enough to the stable point of the map $g_\lambda$, and if n is large, then the trajectory of the avalanche model will stay close to the deterministic sequence $\psi_k$ for a time which is exponentially large in the network size n. For the reader’s convenience we give a short proof of the theorem which generally follows the line of argument in [Reference Aydogmus6] but differs in several details.

Theorem 6.4. Suppose that Assumption 3.1 is satisfied with $\lambda>1.$ For $\delta>0,$ let

(6.12) \begin{eqnarray}\tau_n(\delta)=\inf\bigl\{k\in\mathbb{Z}_+:|x_{n,k}-\psi_k|\geq \delta\bigr\},\end{eqnarray}

where the sequence $\psi_k$ is defined in (6.7).

There exist an interval $(a,b)\subset (0,1)$ including $\zeta_\lambda$ and constants $\gamma >0,$ $\delta_0\in (0,1),$ and $n_0\in\mathbb{N}$ such that if $x_{n,0}\in (a+\delta_0,b-\delta_0),$ and $n>n_0,$ then for any $m\in\mathbb{N}$ and $\delta \in (0,\delta_0)$ we have

\begin{eqnarray*}\mathbb{P}(\tau_n(\delta)>m)\geq \bigl(1-2e^{-\gamma \delta^2 n}\bigr)^m\geq 1-2me^{-\gamma \delta^2 n}.\end{eqnarray*}

Furthermore, the constants a, b, $\delta_0$, and $\gamma$ depend on the sequence of parameters $(p_n)_{n\in\mathbb{N}}$ only through $\lambda$ (this is not necessarily true for $n_0,$ which in general is not exclusively determined by the value of $\lambda$).

Proof. Pick first $b\in (0,1)$ and then $a\in (0,1)$ in such a manner that

\begin{eqnarray*}0<a<\min\{\nu_\lambda,\zeta_\lambda,g_\lambda(b)\}<\max\{\nu_\lambda,\zeta_\lambda\}<b<1.\end{eqnarray*}

Let $I=(a,b)$ and $h=\min\{g_\lambda(a),g_\lambda(b)\}.$ Then

\begin{eqnarray*}a<h,\qquad\nu_\lambda<b,\qquad g(I)\subset (h,\nu_\lambda).\end{eqnarray*}

Thus, if we set

\begin{eqnarray*}\varepsilon=\min\bigl\{b-\zeta_\lambda,h-a,(b-a)/2\bigr\},\end{eqnarray*}

we get

(6.13) \begin{eqnarray}g(a,b)\subset (a+\varepsilon,b-\varepsilon)\qquad \textrm{and}\qquad \zeta_\lambda\in (a+\varepsilon,b-\varepsilon).\end{eqnarray}

The latter assertion is true because $g_\lambda(\zeta_\lambda)=\zeta_\lambda$ and the point $\zeta_\lambda$ belongs to the interval (a,b), which is mapped into $(a+\varepsilon,b-\varepsilon)$ by $g_\lambda.$

Next, observe that for any $x\in(0,1),$

\begin{eqnarray*}g_\lambda^{\prime}(x)=-1+(1+\alpha-\alpha x)e^{-\alpha x}>-1,\end{eqnarray*}

and

\begin{eqnarray*}g_\lambda^{\prime\prime}(x)<0 \quad (\mbox{so that $g_\lambda^{\prime}$ is decreasing}),\qquad g_\lambda^{\prime}(1)=-1+e^{-\alpha}<0.\end{eqnarray*}

Furthermore,

\begin{eqnarray*}g_\lambda^{\prime}(\zeta_\lambda)<1,\end{eqnarray*}

because $g_\lambda^{\prime}(\zeta_\lambda)<0$ when $\zeta_\lambda>\nu_\lambda,$ and when $\zeta_\lambda<\nu_\lambda$ the graph of $g_\lambda$ intersects the line $y=x$ at $\zeta_\lambda$ going upward from left to right, which means the slope is less than one at the point of intersection.

It follows that we can choose $a,b\in(0,1)$ in such a way that (6.13) holds for some $\varepsilon>0,$ and furthermore the following holds:

(6.14) \begin{eqnarray}\textrm{There exists}\ \varrho\in (0,1)\ \textrm{such that}\ |g_\lambda^{\prime}(x)|<\varrho\ \textrm{on}\ (a,1).\end{eqnarray}

From now on assume that the constants $a,b\in (0,1)$ and $\varepsilon>0$ satisfy (6.13) and (6.14). Pick any $\delta\in (0,\varepsilon),$ and assume that for some $n\in\mathbb{N}$ and $k\in\mathbb{Z}_+,$ $X^{(n)}_k$ satisfies the following two conditions:

(6.15) \begin{eqnarray}\begin{array}{l@{\quad}l@{\quad}l}1.&\mbox{}\quad& x_{n,k}\in (a+\delta,b-\delta)\\[3pt]2.&\mbox{}\quad&|x_{n,k}-\psi_k| < \delta.\end{array}\end{eqnarray}

It follows from (1.2) that for a given $X^{(n)}_k<bn,$

(6.16) \begin{eqnarray}\nonumber&&\mathbb{P}\Bigl(\Bigl|x_{n,k+1}-(1-x_{n,k})\Bigl(1-q_n^{X^{(n)}_k}\Bigr)\Bigr|\geq \frac{(1-\varrho)\delta}{2}\,\Bigr|X^{(n)}_k\Bigl)\leq 2e^{-\frac{\delta^2 n(1-x_{n,k})}{2(1-\varrho)^2}}\\&&\qquad\leq 2e^{-\frac{\delta^2 n(1-b)}{2(1-\varrho)^2}},\end{eqnarray}

where in the first step we applied Hoeffding’s inequality for binomial variables. Furthermore, if n is large enough, then under the condition (6.15), we have

\begin{eqnarray*}\Bigl|g_\lambda(x_{n,k})-\bigl(1-x_{n,k}\bigr)\Bigl(1-q_n^{X^{(n)}_k}\Bigr)\Bigr|\leq \frac{(1-\varrho)\delta}{2}\end{eqnarray*}

and

\begin{eqnarray*}|g_\lambda(x_{n,k})-g_\lambda(\psi_k)|\leq \varrho \delta,\end{eqnarray*}

which together imply

\begin{eqnarray*}\Bigl|g_\lambda(\psi_k)-\bigl(1-x_{n,k}\bigr)\Bigl(1-q_n^{X^{(n)}_k}\Bigr)\Bigr|\leq \frac{(1-\varrho)\delta}{2}+\varrho \delta\leq \frac{(1+\varrho)\delta}{2}.\end{eqnarray*}

Combining the last inequality with (6.16), we obtain that

\begin{eqnarray*}\mathbb{P}\Bigl(\bigl|x_{n,k+1}-\psi_{k+1}\bigr|\geq \delta\,\Bigl|\,X^{(n)}_k \Bigl)\leq 2e^{-\frac{\delta^2 n(1-b)}{2(1-\varrho)^2}}\end{eqnarray*}

for any $X_{n,k}$ that satisfies condition (6.15). Taking into account (6.13), we arrive at the following conclusion.

Lemma 6.2. If $X^{(n)}_k$ satisfies the condition (6.15), then (conditionally on $X^{(n)}_k$) $X^{(n)}_{k+1}$ satisfies the same condition with probability larger than

\begin{equation*}2e^{-\frac{\delta^2 n(1-b)}{2(1-\varrho)^2}},\end{equation*}

uniformly on $X^{(n)}_{k+1}.$

Let $A_{n,m,\delta}$ be the event that (6.15) is satisfied for $k=0,1,\ldots,m,$ and recall $\tau_n(\delta)$ from (6.12). By the Markov property, the lemma implies that under the conditions of the theorem, if (6.15) is satisfied for $k=0,$ then

\begin{eqnarray*}\mathbb{P}\bigl(\tau_n(\delta) > m\bigr)\geq \mathbb{P}(A_{n,m,\delta})\geq \Bigl(1-2e^{-\frac{\delta^2 n(1-b)}{2(1-\varrho)^2}}\Bigr)^m\geq 1-2me^{-\frac{\delta^2 n(1-b)}{2(1-\varrho)^2}},\end{eqnarray*}

completing the proof of the theorem.

The counterpart of Theorem 6.2 for $\lambda\leq 1$ is easier to prove, since the behavior of the derivative $g_\lambda^{\prime}$ on [0, 1] is ‘friendlier’ in this case, in that $|g_\lambda^{\prime}(x)|<1$ and $g_\lambda(x)<x$ for all $x\in(0,1).$ For $\varepsilon\in (0,1)$ let

(6.17) \begin{eqnarray}\varsigma_n(\varepsilon)=\inf\bigl\{k\in\mathbb{Z}_+:x_{n,k}<\varepsilon\bigr\}.\end{eqnarray}

Recall $\tau_n(\delta)$ from (6.12). Let $\lfloor x \rfloor =\max\{k\in\mathbb{Z}: k\leq x\}$ denote the integer part of the real number x. We have the following result.

Theorem 6.5. Let Assumption 3.1 hold with $\lambda\leq 1,$ and suppose that $x_{n,0}=\psi_0$ for some $\psi_0\in(0,1)$ and all $n\in\mathbb{N}.$ Pick any $\delta\in (0,1)$ such that

\begin{eqnarray*}\delta<\psi_0-\psi_1,\end{eqnarray*}

where the $\psi_k$ are defined in (6.7). Let $a\in(0,1)$ be the minimal root of the equation $g_\lambda(a)=\delta,$ and define

(6.18) \begin{eqnarray}\varrho=\left\{\begin{array}{l@{\quad}l@{\quad}l}\max\{\lambda, |g_\lambda^{\prime}(\psi_0)|\}&\textrm{if}&\lambda<1,\\ \\[-7pt]\max\{g^{\prime}_\lambda(a), |g_\lambda^{\prime}(\psi_0)|\}&\textrm{if}&\lambda=1.\end{array}\right.\end{eqnarray}

Then the following hold:

1. (i) Let $m_0=\bigl\lfloor \log \frac{\psi_0}{\delta}\bigr\rfloor +1.$ Then for any $n\in\mathbb{N},$

   \begin{eqnarray*}\mathbb{P}(\varsigma_n(\delta)>m_0)\geq \left(1-2e^{-\frac{\delta^2 n(1-\psi_0)}{2(1-\varrho)^2}}\right)^{m_0}\geq 1-2m_0\exp\Bigl\{-\frac{\delta^2 n(1-\psi_0)}{2(1-\varrho)^2}\Bigr\}.\end{eqnarray*}

2. (ii) If $\lambda\in (0,1),$ then for any $n,m\in\mathbb{N},$

   \begin{eqnarray*}\mathbb{P}(\tau_n(\delta)>m)\geq \left(1-2e^{-\frac{\delta^2 n(1-\psi_0)}{2(1-\varrho)^2}}\right)^m\geq 1-2m\exp\Bigl\{-\frac{\delta^2 n(1-\psi_0)}{2(1-\varrho)^2}\Bigr\}.\end{eqnarray*}

Proof. Observe that (6.14) holds with $\varrho$ introduced in (6.18). Furthermore, $\psi_k$ is a decreasing sequence since $g_\lambda(x)<x$ for $x\in (0,1)$ when $\lambda\leq 1.$ The rest of the proof is similar to that of Theorem 6.4; specifically, the induction argument based on (6.15) and Lemma 6.2 carries over. Note that we can replace $1-b$ by $1-\psi_0=1-\varepsilon_0$ in the conclusions of the lemma, because the sequence $\psi_k$ is monotone decreasing under the conditions of the theorem.