Structure of the particle population for a branching random walk with a critical reproduction law

TL;DR

This paper investigates the structure of particle populations in a symmetric branching random walk on a lattice, revealing how subpopulation behaviors influence convergence and clustering in low dimensions.

Contribution

It provides new insights into the subpopulation structure and explains the effects of vanishing subpopulations on steady state convergence and cluster formation.

Findings

Vanishing subpopulations do not hinder convergence to steady state.

Clusterization occurs in dimensions 1 and 2 due to subpopulation dynamics.

The structure of subpopulations explains the behavior of the entire process.

Abstract

We consider a continuous-time symmetric branching random walk on the $d$-dimensional lattice, $d\ge 1$, and assume that at the initial moment there is one particle at every lattice point. Moreover, we assume that the underlying random walk has a finite variance of jumps and the reproduction law is described by a critical Bienamye-Galton-Watson process at every lattice point. We study the structure of the particle subpopulation generated by the initial particle situated at a lattice point $x$. We answer why vanishing of the majority of subpopulations does not affect the convergence to the steady state and leads to clusterization for lattice dimensions $d=1$ and $d=2$.

Equations

$$\mathcal{L}\psi (x)=\kappa \sum _{z\ne 0}[\psi (x+z)-\psi (x)]a(z).$$

$$\sum _{z\in {\mathbb{Z}}^d}{e}^{(\lambda ,z)}a(z)<\infty .$$

$$\frac{\partial p}{\partial t}=\mathcal{L}p,p(0,x,y)=\delta (y-x).$$

$$\hat{a}(k)=\sum _{z\ne 0}{e}^{i(k,x)}a(z)=\sum _{z\ne 0}\cos (kz)a(z).$$

$$p(t,x,y)=\frac{1}{(2\pi {)}^d}\underset{T^d}{\int }{e}^{-\kappa t(1-\hat{a}(k)+ik(x-y))}dk,T^d=[-\pi ,\pi {]}^d.$$

$$p(t,0,z)\le p(t,0,0)=\frac{1}{(2\pi {)}^d}\underset{T^d}{\int }{e}^{-\kappa t(1-\hat{a}(k))}dk.$$

$$p(t,z)\sim \frac{e^{-\frac{|z{|}^2}{2t}}}{(2\pi t{)}^{d/2}},t\to \infty .$$

$$n(t,y)=\sum _{x\in {\mathbb{Z}}^d}n(t,x,y),$$

$$K_t(x_1,...,x_m)=\mathsf{E}n(t,x_1)...n(t,x_m),$$

$$\frac{\partial u_z(t,x,y)}{\partial t}$$

$$u_z(0,x,y)$$

$$\frac{\partial u_z(t,x,y)}{\partial t}$$

$$u_z(0,x;y)$$

$$m_{k}(t,x,y):=\mathsf{E}\prod_{i=0}^{k-1}(n(t,x,y)-i)=\frac{\partial^{k}u_{z}(t,x,y)}{\partial z^{k}}\bigg{|}_{z=1}.$$

$$\frac{\partial m_1(t,x,y)}{\partial t}$$

$$m_1(0,x,y)$$

$$m_1(t,x,y)=p(t,x,y)$$

$$\mathsf{E}n(t,y)=\sum _{x\in {\mathbb{Z}}^d}\mathsf{E}n(t,x,y)=\sum _{x\in {\mathbb{Z}}^d}{m}_1(t,x,y)=\sum _{x\in {\mathbb{Z}}^d}p(t,x,y)\equiv 1$$

$$\frac{\partial m_2(t,x,y)}{\partial t}$$

$$m_2(0,x,y)$$

$$\mathcal{F}(\mathcal{L}(m_1(t,x,y)))$$

$$\hat{\mathcal{L}}(k)$$

$$\frac{\partial {\hat{m}}_1(t,x,y)}{\partial t}$$

$${\hat{m}}_1(0,k,y)$$

$${\hat{m}}_1(t,k,y)=e^{t\hat{\mathcal{L}}(k)}{e}^{i(k,y)},$$

$$\frac{\partial {\hat{m}}_2(t,k,y)}{\partial t}$$

$${\hat{m}}_2(0,k,y)$$

$${\hat{m}}_2(t,k,y)e^{-t\hat{\mathcal{L}}(k)}=2\beta {\int }_0^t{\hat{m}}_1(t,k,y)\ast {\hat{m}}_1(t,k,y)e^{-s\hat{\mathcal{L}}(k)}ds,$$

$$m_k(t,x,\Gamma ):=\sum _{y\in \Gamma }{m}_k(t,x,y)$$

$$\sum _{x\in \mathbb{Z}}{m}_2(t,x,\Gamma )$$

$$=2\beta {\int }_{T^d}d\theta {\int }_0^t\sum _{y_1\in \Gamma }\sum _{y_2\in \Gamma }\cos ((\theta ,y_2-y_1))e^{2s\hat{\mathcal{L}}(\theta )}ds$$

$$=\beta {\int }_{T^d}\frac{(e^{2t\hat{\mathcal{L}}(\theta )}-1)(\sum _{y_1\in \Gamma }\sum _{y_2\in \Gamma }\cos ((\theta ,y_2-y_1)))}{\hat{\mathcal{L}}(\theta )}d\theta .$$

$$G_0(0,0)$$

$$=\frac{1}{(2\pi {)}^d}\underset{T^d}{\int }\frac{dk}{\kappa (1-\hat{a}(k))}=\infty ,$$

$$m_2(t,\Gamma )\sim {\hat{m}}_2(t,0,\Gamma )\to \beta {\int }_{T^d}\frac{\sum _{y_1\in \Gamma }\sum _{y_2\in \Gamma }\cos ((\theta ,y_2-y_1))d\theta }{-\hat{\mathcal{L}}(\theta )}\to \infty ,$$

Full text

Keywords: Branching random walk; critical branching process; limit theorems; population dynamics.

MSC 2010: 60J80, 60J35, 60F10

Structure of the particle population for a branching random walk with a critical reproduction law

Daria Balashova, Stanislav Molchanov, Elena Yarovaya

1 Introduction

Consider a lattice population model, i.e. a random field $n(t,\cdot)$ of particles on $\mathbb{Z}^{d}$, $d\geq 1$, where $n(t,y)$ is the number of particles at the point $y\in\mathbb{Z}^{d}$ at the time moment $t\geq 0$. Let $n(0,y)\equiv 1$ for every $y\in\mathbb{Z}^{d}$. The spatio-temporal evolution of the field includes the migration and birth-death processes. As usual, we exclude interaction between particles.

The migration until the first transformation (death or particle splitting) is described by a random walk of particles governed by the generator

[TABLE]

Here $\kappa>0$ is the diffusion coefficient and $a(z)\geq 0$, $\sum\limits_{z\neq 0}a(z)=-a(0)=1$ is the distribution of the random walk jumps. We assume symmetry $a(z)=a(-z)$ and irreducibility of the random walk, see [5]. Finally, we assume that, for every $\lambda\in\mathbb{R}^{d}$,

[TABLE]

It means that the tails of $a(z)$ are superexponentially light. Note that the large deviation theorems from [5] based on the last condition.

The transition probabilities of the random walk $X(t),t\geq 0$, with the generator $\mathcal{L}$, i.e. $p(t,x,y)=P_{x}\{X(t)=y\}$, satisfy, see [1], the Kolmogorov backward equations

[TABLE]

Consider the Fourier transform

[TABLE]

Due to the homogeneity of the random walk over the space, we have $p(t,x,y)=p(t,0,y-x)=p(t,0,x-y)$. By applying the Fourier transform to equation (2) we get

[TABLE]

Note that $\hat{a}(k)$ is a real function and as a result we have the following important inequality:

[TABLE]

In what follows we use the notation: $p(t,z):=p(t,0,z)$ and $p(t,y-x):=p(t,x,y)$.

We assume that the intensity of jumps $\kappa$ in (1) is equal to one, then the local Central Limit Theorem (CLT) from [5] implies for $|z|<C\sqrt{t}$ that

[TABLE]

Each particle of the population dies in the interval ($t,t+dt$) with probability $\mu dt$, where $\mu$ is the mortality rate, or splits into two identical particles with probability $\beta dt$. We call $\beta$ the birth rate. Further we consider only the critical branching process: $\beta=\mu$ at every lattice point.

The total population $n(t,y)$ at a point $y\in\mathbb{Z}^{d}$ is the sum of independent subpopulations:

[TABLE]

where $n(0,x,y)=\delta(y-x)$, and $n(0,y)\equiv 1$.

For the moment analysis of the field $n(t,y)$ one can use either the forward Kolmogorov equations for the correlation functions

[TABLE]

or the backward approach for $n(t,x,y)$ with the following summation over $x$ of the cumulants of subpopulations. Note that the first approach is simpler for continuous models with $\mathbb{R}^{d}$ instead of $\mathbb{Z}^{d}$, see [2]. Their method requires one particle in each site and is not applicable to the lattice model. In the latter case the backward approach was developed in [3].

The papers [2], [3] contain the following results: a branching random walk for $\mu=\beta$ on $\mathbb{Z}^{d}$ (or $\mathbb{R}^{d}$) converges to a steady state (statistical equilibrium) if the underlying process $X(t)$ with the generator $\mathcal{L}$ is transient.

Remark 1**.**

Theorems 3.1 and 3.2 in [2] contain some additional assumptions about the function $a(x)$, $x\in\mathbb{R}^{d}$. In the lattice case we have no additional assumptions on $X(t)$, except the transience.

For the generating function of subpopulation $u_{z}(t,x,y)=\mathsf{E}z^{n(t,x,y)}$ we have the parabolic problem with quadratic non-linearity (see (16) in [3]):

[TABLE]

From this last equation in the critical case $\mu=\beta$ we get

[TABLE]

We define the factorial moments as

[TABLE]

Then, differentiating the equation (6) with respect to $z$ and substituting $z=1$, we obtain

[TABLE]

Note that the Cauchy problems for $p(t,x,y)$ (2) and $m_{1}(t,x,y)$ are the same, so we get

[TABLE]

and as a result

[TABLE]

(conservation law for the density of the total population).

The analysis of the second and higher moments ($k\geq 3$) for the total population and subpopulations is more complex, see details in [3].

For the second moment $m_{2}(t,x,y)=\mathsf{E}[n(t,x,y)(n(t,x,y)-1)]$ we have the following results:

[TABLE]

We consider the Fourier representations

[TABLE]

Thus, one can obtain

[TABLE]

its solution takes the form

[TABLE]

and for

[TABLE]

we obtain

[TABLE]

where “$\ast$” denotes the convolution operator.

Denote

[TABLE]

and note that

[TABLE]

If the random walk generated by the operator $\mathcal{L}$ is recurrent, i.e.

[TABLE]

then for $t\to\infty$

[TABLE]

At the physical level this is a manifestation of the high irregularity of the field $n(t,y)$. For large $t$ such a field (in which the density is meant to be 1!) consists of big islands separated by large distances (the phenomenon which has different names: clusterization, clumping, intermittency, see [6, 7]). Our goal is to explain this phenomenon.

Note that the recurrence of $X(t)$ with a finite variance of jumps can appear only in dimensions $d=1$ and $d=2$, and implies that

[TABLE]

If, however, the process $X(t)$ is transient, i.e

[TABLE]

then field $n(t,y)$ converges in law to a steady state (statistical equilibrium), see [3].

On the other hand, if we consider the sums

[TABLE]

i.e. the total number of particles in the subpopulation, then they are independent Galton-Watson processes and the generating function

[TABLE]

is the solution of the differential equation

[TABLE]

Proposition 1**.**

[TABLE]

Proof.

By integrating (7) we obtain

[TABLE]

Therefore,

[TABLE]

Then

[TABLE]

∎

Theorem 1**.**

If $t\to\infty$ and $s>0$ then

[TABLE]

Proof.

[TABLE]

∎

We have the following situation: the majority of the subpopulations $n_{x}(t)$ are vanishing on the large time interval $[0,t]$ but for the remaining subpopulations, proportion of which is $\sim(\beta t)^{-1}$, the number of particles will have the order $O(t)$ at least at the level of the first moment. What is the structure of such large conditional subpopulations? We have to answer this question if we want to understand why vanishing of majority of the subpopulations does not affect the convergence to the steady state in the recurrent case $G_{0}(0,0)=\infty$ and leads to clusterization.

2 The structure of the subpopulation $n(t,x,y)$, $y\in\mathbb{Z}^{d}$, for fixed $x$, large $t$, and $n_{x}(t)>0$

The joint generating function of the processes $n_{x}(t)$ and $n(t,x,y)$, i.e.

[TABLE]

is the solution of the same equation (5) but with the different critical condition:

[TABLE]

We cannot solve the non-linear equation (12) but using the differentiation over $z_{1}$ and substitution $z_{1}=1$ one can calculate or estimate the conditional moments $\mathsf{E}[n(t,x,y)|n_{x}(t)>0]$, $\mathsf{E}[n(t,x,y)(n(t,x,y)-1)|n_{x}(t)>0]$ etc.

Put

[TABLE]

Denote $\gamma(t,z):=2\beta(u_{z,z_{1}=1}(t,x,y)-1)$, then

[TABLE]

Note that due to (9) we have

[TABLE]

Lemma 1**.**

Consider the parabolic problem with the potential, depending only on the time $t$:

[TABLE]

Then

[TABLE]

Proof.

Consider

[TABLE]

then

[TABLE]

Taking a Fourier transform gives:

[TABLE]

Note that

[TABLE]

then

[TABLE]

Therefore, we obtain

[TABLE]

∎

Lemma 2**.**

[TABLE]

Proof.

Let $X,Y$ be two integer-valued nonnegative random variables then

[TABLE]

Then

[TABLE]

Let us return to the equation (13). Note that

[TABLE]

then

[TABLE]

We used here the expansion of \bigg{(}1-\frac{\beta tz}{\beta t+1}\bigg{)}^{-2} into the Taylor series. Since

[TABLE]

we have finally that

[TABLE]

∎

Theorem 2**.**

For $t\to\infty$ we have

[TABLE]

Proof.

[TABLE]

Then using (1) we have

[TABLE]

∎

Then for $|x-y|=O(\sqrt{t})$

[TABLE]

The formula (14) is not difficult to understand without calculations. Each particle, among $n_{x}(t)=k$ particles, performs on $[0,t]$ a random walk with the generator $\mathcal{L}$ (which is the union of the pieces of paths between the successive transformations). These paths are highly dependent, but for calculating the first moment this dependence is irrelevant.

For the second (factorial) moment we can use the same approach. Differentiating (12) two times over $z_{1}$ and substituting $z_{1}=1$ we get for

[TABLE]

the equation

[TABLE]

Here

[TABLE]

is the same space independent potential as in $\tilde{m}_{1}(t,x,y,z)$, see above.

Lemma 3**.**

Consider the non-homogeneous parabolic problem

[TABLE]

Then the substitution $u=e^{\int\limits_{0}^{t}\gamma(s)\,ds}v(t,x)$ leads to the equation

[TABLE]

and its solution is given by Duhamel’s formula

[TABLE]

Proof.

Taking a Fourier transform gives:

[TABLE]

Note that

[TABLE]

then

[TABLE]

∎

Combining together formulas (15), (16) and expressions

[TABLE]

we can find the conditional second moment

[TABLE]

Let us now represent the integrand by the Taylor series. Since

[TABLE]

we will get

[TABLE]

i.e. (setting $m=n+2$) we will get finally

[TABLE]

Due to lemma 2 we have the following formula for the conditional second moment

[TABLE]

Like in the case of the first moment the most important values of $m$ must be of order $t$. To simplify the calculations we will put $m=t,\beta=1$. Transition to the more general case $m=ct,\beta\neq 1$ is simple.

We start from the rough estimation of the second moment $M_{2}(t,t,\cdot)$. The const $c$ in the calculations below will be “universal” , i.e. its meaning can be different in the neighbouring formulas. We start from several simple inequalities.

[TABLE]

Then formula (17) together with four inequalities above for all $x,y$ gives

[TABLE]

It gives

[TABLE]

Assume $a(z)$ is bounded with the following asymptotic expansion at infinity:

[TABLE]

where $\alpha>2$, $\alpha_{0}=\alpha<\alpha_{1}<\ldots<\alpha_{N}$, $\dot{z}=z/|z|$, and $c_{0}(\dot{z})=c_{0}(\dot{-z})$ is a positive continuous function on $\mathcal{S}^{d-1}$, $l=1$ if $\alpha>[\alpha]$, $l=2$ if $\alpha=[\alpha]$, $c_{0}(\dot{z})>0$, and $c_{j}(\dot{z})$ are sufficiently smooth.

Theorem 3**.**

Let condition (19) be satisfied, for $|x-y|<C\sqrt{t}$ and $t\to\infty$

[TABLE]

Proof.

From ([4], Theorem 1.1) there are constants $A$, $\epsilon>0$, such that the solution $p(t,z)$ has the following asymptotic behavior as $|z|^{2}\geq At$, $t\geq 0$, $x\to\infty$:

[TABLE]

and from the local CLT

[TABLE]

then

[TABLE]

[TABLE]

Due to (17) we obtain

[TABLE]

Denote $q=\frac{s(\beta t+1)}{t(\beta s+1)}$, then for $m=t$ and $\beta=1$

[TABLE]

then

[TABLE]

For $d=1$

[TABLE]

for $d=2$

[TABLE]

for $d\geq 3$ and constants $q_{i}(i=0,\ldots,d-1$) and $r_{1},r_{2}$ the integrand $M_{2}(t,x,y)$ has following form: for $d=2n+1$

[TABLE]

and for $d=2n$

[TABLE]

Therefore, one can combine this inequalities with (18) to obtain (20). ∎

3 Conclusion

Assuming that $n(0,x)\equiv 1$ for the large moment $t$ the majority of the subpopulation $n(t,x,\cdot)$ will degenerate, see (8) and (10). The initial points of all such (non-degenerated) subpopulations form the Bernoulli point field with parameter $p=\frac{1}{1+\beta t}$.

In a neighborhood of each such point the subpopulation forms the island $n(t,x,\cdot)$ containing the random number $n_{x_{i}}(t)$ of the particles with the mean value $\beta t+1$. Then Theorem 3.1 from [4] gives

[TABLE]

i.e. the diameter of the subpopulation is $O(\sqrt{t})$.

Consider several cases: $d=1$, $d=2$ and $d\geq 3$.

Case $d=1$.

Here the distances between the points $x_{i,t}$ have the geometric distribution with the mean value $\beta t$ and the subpopulations have typical size $\sqrt{t}$ (due to local CLT).

The total population demonstrates the high level of intermittency: large clusters, that is clusters of the diameter $O(\sqrt{t})$, are separated by the empty intervals of the length $O(t)$. Due to (8) each of such clusters contains about $t$ particles.

Case $d=2$.

Here the typical size of the subpopulation and the typical distance between points $x_{i,t}$ both have order $\sqrt{t}$. However, the population still has fairly large gaps. Let us estimate these gaps.

Consider on the lattice $\mathbb{Z}^{2}$ the square $[-L,L]\times[-L,L]=Q_{L}$ with large enough $L$, containing $(2L+1)^{2}$ point. For $t=0$ any point $x\in Q_{L}$ contains one particle: at the moment $t\gg 1$ only $\frac{(2L+1)^{2}}{\beta t}=N_{t}$ subpopulations started from the points $x_{i,t}$ will survive.

Let us divide $Q_{L}$ into cells of the size ($\sqrt{t}a(t)\times\sqrt{t}a(t)$) where the function $a(t)$ will be selected later. Number of such cells is $\frac{(2L+1)^{2}}{ta^{2}(t)}=N_{1}(t)$. The probability that one cell does not contain points $x_{i,t}$ equals to

[TABLE]

Then the mean number of empty cells in $Q_{L}$ will be

[TABLE]

Put $a^{2}(t)=\beta\ln t$. Then $a(t)\sim c\sqrt{\ln t}$ for large $t$. Hence for $L=t\sqrt{\ln t}$ we obtain $\mu_{t}\sim\mathrm{const}$, i.e. empty cells have density $O(\frac{1}{t\sqrt{\ln t}})$, which means that in each cube of the size $t\sqrt{\ln t}\times t\sqrt{\ln t}$ there is a gap of the diameter $O(t\sqrt{\ln t})$. This is a very weak intermittency.

Case $d\geq 3$.

Here the population $n(t,x)$ is highly uniform.

4 Simulation

We call the state of a BRW system the set of pairs $\{x,t\}$, each of which corresponds to a particle located at the point $x\in\mathbb{Z}^{d}$, that first appeared at the point $x$ at the time $t$. By evolution we mean a jump to another point, splitting or death of a particle. In the process of modeling the transition from one state of the BRW system to another will be carried out by excluding one pair from the set of states and adding to the seat of states of one or several pairs corresponding to the result of the simulated particle evolution.

Initialization. First, we set the characteristics of the simulated BRW: choose the dimension $d$ of the integer lattice, the functions defining the distribution of the jumps matrix $A$ and the infinitesimal generating functions conditioned by the relation $\beta=\mu$, as well as the execution time $T$ and a sufficiently large part of the space in which we will consider the process. At the initial moment of time, the state of the system is determined by the presence of a single particle at each point of a given space.

Step of algorithm. On input we have pairs $\{x,t\}$ waiting to be processed. We select an arbitrary pair $\{x,t\}$ (due to the independence of the particles) and model the evolution of the corresponding particle: the exponential time $dt$ of staying it in $x$, then the jump or birth/death event. In the case of a jump, the transition state is simulated according to the matrix $A$, in the next state the current pair $\{x,t\}$ disappears and appears new pair $\{x^{\prime},t+dt\}$. In the case of death, the current pair $\{x,t\}$ disappears, and in the case of birth (dividing into two offsprings) the current pair $\{x,t\}$ disappears and two new pairs $\{x,t+dt\}$ are added to the set.

Stop condition. The algorithm terminates when all the values of $t$ of all pairs in the system exceed the specified time $T$, or when the set of states becomes empty (the process has degenerated). Finally, we obtain the coordinates of all the points in the system remaining at time $T$ and plot them on the graph.

At a point $x\in\mathbb{Z}^{d}$ the particle can perform a symmetric random walk or branching, which is a Galton-Watson process with infinitesimal generating function $f(u)=\beta-2\beta u+\beta u^{2}$. In this case, the particle’s behavior is described as follows: the particle spends at a point $x$ a random time, exponentially distributed with the parameter $1+2\beta$, and then goes to the point $x^{\prime}$ with probability $a(x,x^{\prime})/(1+2\beta)$ or die with probability $\beta/(1+2\beta)$ or divided into two offsprings with probability $\beta/(1+2\beta)$.

To simulate “randomness” in implementing the proposed algorithm, we use the following random number generators:

[TABLE]

In the case of a walk, the transition probabilities from $x$ to $x^{\prime}$ are given by

[TABLE]

and in the case of splitting/death

[TABLE]

For $\mathbb{Z}^{1}$ the obtained clusters are plotted in Fig. 1 and for $\mathbb{Z}^{2}$ they are plotted in Fig. 2.

Acknowledgements. D. Balashova and E. Yarovaya were supported by the Russian Foundation for Basic Research (RFBR), project No. 17-01-00468. S. Molchanov was supported by the Russian Science Foundation (RSF), project No. 17-11-01098.