Population Processes with Immigration

TL;DR

This paper provides a comprehensive analysis of Galton-Watson processes with immigration, including models in random environments and branching random walks, establishing key limit results for correlation functions.

Contribution

It offers the first complete analysis of Galton-Watson models with immigration across various environments and extends to branching random walks with new limit theorems.

Findings

Existence of limits for the first two correlation functions in branching random walks with immigration

Complete characterization of Galton-Watson models with immigration in different environments

Analysis includes stationary and non-stationary processes

Abstract

The paper contains the complete analysis of the Galton-Watson models with immigration, including the processes in the random environment, stationary or non-stationary ones. We also study the branching random walk on $Z^d$ with immigration and prove the existence of the limits for the first two correlation functions.

Figures

Equations

$${\phi }_t(\lambda )=E{e}^{-\lambda n(t,x)}=\sum _{j=0}^{\infty }P\{n(t,x)=j\}{e}^{-\lambda j}.$$

$$m_j(t,x)=(-1{)}^j\frac{{\partial }^j\phi }{\partial {\lambda }^j}{|}_{\lambda =0}.$$

$$n(t+dt,x)=n(t,x)+{\xi }_{dt}(t,x)$$

$$\xi_{dt}(t,x)=\left\{\begin{array}[]{l}+1\hskip 14.22636pt\beta n(t,x)dt+kdt\\ -1\hskip 14.22636pt\mu n(t,x)dt\\ 0\hskip 22.76228pt1-((\beta+\mu)n(t,x)+k)dt\\ \end{array}\right.$$

$$\frac{\partial {\phi }_t(\lambda )}{\partial t}={\phi }_t(\lambda )\left(\mu n(t)e^{\lambda }-((\beta +\mu )n(t)+k)+(\beta n(t)+k)e^{-\lambda }\right)$$

$${\phi }_0(\lambda )=e^{-\lambda }$$

$$\frac{\partial {\phi }_t(\lambda {)}^{(j)}}{\partial t}={\phi }_t(\lambda {)}^{(j)}\left(\mu n(t)e^{\lambda }-((\beta +\mu )n(t)+k)+(\beta n(t)+k)e^{-\lambda }\right)+$$

$$+\sum _{i=1}^j\left(\genfrac{}{}{0ex}{}{j}{i}\right){\phi }_t(\lambda {)}^{(j-i)}\left(\mu n(t)e^{\lambda }+(-1{)}^i(\beta n(t)+k)e^{-\lambda }\right)$$

$${\phi }_0(\lambda {)}^{(j)}=(-1{)}^j{e}^{-\lambda }$$

$$\frac{d{m}_j(t)}{dt}=\sum _{i=1}^j\left(\genfrac{}{}{0ex}{}{j}{i}\right)\left((\beta +(-1{)}^i\mu )m_{j-i+1}+m_{j-i}\right)$$

$$=j(\beta -\mu )m_j+s_j$$

$$m_j(0)=1$$

$$\frac{d{m}_1(t)}{dt}=(\beta -\mu )m_1(t)+k$$

$$m_1(0)=1$$

$$\frac{d{m}_2(t)}{dt}=2(\beta -\mu )m_2(t)+(\beta +\mu +2k)m_1(t)+k$$

$$m_2(0)=1$$

$$m_1(t)=\frac{k}{\mu -\beta }+(1-\frac{k}{\mu -\beta })e^{-(\mu -\beta )t}$$

$$\begin{array}{rl}{m}_2(t)=& \frac{k(k+\mu )}{(\mu -\beta {)}^2}+\frac{{\mu }^2-2k^2-{\beta }^2+k\mu -3k\beta }{(\mu -\beta {)}^2}{e}^{-(\mu -\beta )t}+\\ & +\frac{k^2+2{\beta }^2+3k\beta -2\mu \beta -2k\mu }{(\mu -\beta {)}^2}{e}^{-2(\mu -\beta )t}\end{array}$$

$$m_1(t)\underset{t\to \infty }{\to }\frac{k}{\mu -\beta }$$

$$m_2(t)\underset{t\to \infty }{\to }\frac{k(k+\mu )}{(\mu -\beta {)}^2}$$

$$\mathrm{Var}(n(t))=m_2(t)-m_1^2(t)\underset{t\to \infty }{\to }\frac{\mu k}{(\mu -\beta {)}^2}.$$

$$\mathcal{L}\psi (n)$$

$$\mathcal{L}\psi (0)$$

$$\begin{array}[]{ll}\psi_{1}(n)&\equiv 1\\ \psi_{2}(n)&=\left\{\begin{array}[]{ll}0&n=0\\ 1&n=1\\ 1+\frac{\mu_{1}}{\lambda_{1}}+\frac{\mu_{1}\mu_{2}}{\lambda_{1}\lambda_{2}}+\cdots+\frac{\mu_{1}\mu_{2}\cdots\mu_{n-1}}{\lambda_{1}\lambda_{2}\cdots\lambda_{n-1}}&n\geqslant 2\\ \end{array}\right.\\ \end{array}$$

$$\pi (1)$$

$$\pi (2)$$

$$\cdots$$

$$\pi (n)$$

$$x_n=\frac{{\lambda }_0\cdots {\lambda }_{n-1}}{\text{textbullet}{\mu }_1\cdots {\mu }_n}=\frac{k(k+\beta )\cdots (k+(\mu -1))\beta }{\mu (2\mu )\cdots (n\mu )}.$$

$$\pi (n)=\frac{1}{\tilde{S}}\frac{{\lambda }_0\cdots {\lambda }_{n-1}}{\text{textbullet}{\mu }_1\cdots {\mu }_n},$$

$$\tilde{S}=1+\frac{k}{\mu }+\frac{k(\beta +k)}{\mu (2\mu )}+\cdots +\frac{k(k+\beta )\cdots (\beta (n-1)+k)}{\mu (2\mu )\cdots (n\mu )}+\cdots .$$

$$\pi (n_0+l)\sim \frac{e^{-\frac{l^2}{2{\sigma }^2}}}{\sqrt{2\pi {\sigma }^2}}ask\to \infty$$

$$\pi (n)$$

$$=\frac{1}{\tilde{S}}{\left(\frac{\beta }{\mu }\right)}^n\frac{\frac{k}{\beta }(\frac{k}{\beta }+1)\cdots (\frac{k}{\beta }+n-1)}{n!}$$

$$=\frac{1}{\tilde{S}}{\left(\frac{\beta }{\mu }\right)}^n\frac{\Gamma (\frac{k}{\beta }+n)}{\Gamma (\frac{k}{\beta })n!}.$$

$$\tilde{S}={\left(1-\frac{\beta }{\mu }\right)}^{-\frac{k}{\beta }}.$$

$$a_n={\left(\frac{\beta }{\mu }\right)}^n\frac{\Gamma (\frac{k}{\beta }+n)}{\Gamma (\frac{k}{\beta })n!}.$$

Full text

1 Introduction

A problem with many single population models of population dynamics involving processes of birth, death, and migration is that the populations do not attain steady states or do so only under critical conditions. One solution is to allow immigration, which can stabilize the population when the birth rate is less than the mortality rate.

Here, we present analysis of several models that incorporate immigration. The first two are spatial Galton-Watson processes, the first with no migration and the second with finite Markov chain spatial dynamics (see section 2 and 3 respectively). The third model allows migration on $\mathbb{Z}^{d}$ (see section 4). The remaining models all involve random environments in some way (see section 5). Two are again Galton-Watson processes, the first with a random environment based on population size and the second with a random environment given by a Markov chain. The last two models have birth, death, immigration, and migration in a random environment allowing in some way non-stationarity in both space and time. We study in this paper only first and second moments. We will return to the complete analysis of the models with immigration in another publication. It will include a theorem about the existence of steady states and an analysis of the stability of these states.

2 Spatial Galton-Watson process with immigration. No migration and no random environment.

2.1 Moments

Assume that at each site for each particle we have birth of one new particle with rate $\beta$ and death of the particle with rate $\mu$. Also assume that regardless of the number of particles at the site we have immigration of one new particle with rate $k$ (this is a simplified version of the process in [14]). Assume that $\beta<\mu$, for otherwise the population will grow exponentially. Assume we start with one particle at each site. In continuous time, for a given site $x$, $x\in Z^{d}$, we can obtain all moments recursively by means of the Laplace transform with respect to $n(t,x)$, where $n(t,x)$ is the population size at time $t$ at $x$

[TABLE]

Specifically, for the $j$th moment, $m_{j}$

[TABLE]

A partial differential equation for $\varphi_{t}(\lambda)$ can be derived using the forward Kolmogorov equations

[TABLE]

where the r.v.​ $\xi$ is defined

[TABLE]

In other words, our site ($x$) in a small time interval ($dt$) can gain a new particle at rate $\beta$ for every particle at the site or through immigration with rate $k$; it can lose a particle at rate $\mu$ for every particle at the site; or no change at all can happen. Because our model is homogeneous in space, we can write $n(t)$ for $n(t,x)$. This leads to the general differential equation

[TABLE]

from which we can calculate the recursive set of differential equations

[TABLE]

Applying 2.1 we obtain a set of recursive differential equations for the moments

[TABLE]

where $s_{j}$ denotes a linear expression involving lower order moments and where we define $m_{0}=1$. For example, the differential equations for the first and second moments are

[TABLE]

and

[TABLE]

These have the solutions:

[TABLE]

and

[TABLE]

Again, given that we have assumed that $\mu>\beta$, in other words, the birth rate is not high enough to maintain the population size, as $t\to\infty$

[TABLE]

[TABLE]

and

[TABLE]

Moreover, it is clear from Eq. 2.1 that all the moments are finite.

In other words, the population size will approach a finite limit, which can be regulated by controlling the immigration rate $k$, and this population size will be stable, as indicated by the fact that the limiting variance is finite. Without immigration, i.e., if $k=0$, the population size will decay exponentially. Another possibility, because all sites are independent and there are no spatial dynamics, is for there to be immigration at some sites, which therefore reach stable population levels, and not at others, where the population thus decreases exponentially. Of course, if the birth rate exceeds the death rate, $\beta>\mu$, $m_{1}(t)$ increases exponentially and immigration has negligible effect, as shown by the solution for $m_{1}(t)$.

2.2 Local CLT

Setting $\lambda_{n}=n\beta+k$, $\mu_{n}=n\mu$, we see that the model given by Eqs. 2.2 and 2.3 is a particular case of the general random walk on $Z_{+}^{1}=\{{0,1,2,\cdots}\}$ with generator

[TABLE]

The theory of such chains has interesting connections to the theory of orthogonal polynomials, the moments problem, and related topics (see [5]). We recall several facts of this theory.

1. a.

Equation $\mathcal{L}\psi=0,x\geqslant 1$, (i.e., the equation for harmonic functions) has two linearly independent solutions:

[TABLE] 2. b.

Denoting the adjoint of $\mathcal{L}$ by $\mathcal{L}^{*}$, equation $\mathcal{L}^{*}\pi=0$ (i.e., the equation for the stationary distribution, which can be infinite) has the positive solution

[TABLE]

This random walk is ergodic (i.e., $n(t)$ converges to a statistical equilibrium, a steady state) if and only if the series $1+\frac{\lambda_{0}}{\mu_{1}}\cdots+\frac{\lambda_{0}\lambda_{1}}{\mu_{1}\mu_{2}}+\cdots+\frac{\lambda_{0}\lambda_{1}\cdots\lambda_{n-1}}{\mu_{1}\mu_{2}\cdots\mu_{n}}$ converges. In our case,

[TABLE]

If $\beta>\mu$, then, for $n>n_{0}$, for some fixed $\varepsilon>0$, $\frac{k+(n-1)\beta}{n\mu}>1+\varepsilon$, that is, $x_{n}\geq C^{n}$, for $C>1$ and $n\geq n_{1}(\varepsilon)$, and so $\sum x_{n}=\infty$. In contrast, if $\beta<\mu$, then, for some $0<\varepsilon<1$, $\frac{k+(n-1)\beta}{n\mu}<1-\varepsilon$, and $x_{n}\leq q^{n}$, for $0<q<1$ and $n>n_{1}(\varepsilon)$; thus, $\sum x_{n}<\infty$. In this ergodic case, the invariant distribution of the random walk $n(t)$ is given by the formula

[TABLE]

where

[TABLE]

Theorem 2.1** (Local Central Limit theorem).**

Let $\beta<\mu$. If $l=O(k^{2/3})$, then, for the invariant distribution $\pi(n)$

[TABLE]

where $\sigma^{2}=\frac{\mu k}{(\mu-\beta)^{2}}$, $n_{0}\sim\frac{k}{\mu-\beta}$.

Proof.

[TABLE]

We see $\tilde{S}$ is a degenerate hypergeometric series, thus

[TABLE]

Set

[TABLE]

Then, $\pi(n)=\displaystyle\frac{a_{n}}{\tilde{S}}$. We have

[TABLE]

and because $a_{n_{0}}\sim\frac{k}{\mu-\beta}$

[TABLE]

We consider $\sum\limits_{i=1}^{l}\ln(1+\frac{\beta(i-1)(\mu-\beta)}{\mu k})=\int_{1}^{l}\ln(1+\frac{(x-1)(\mu-\beta)\beta}{\mu k})dx+O(\ln(1+\frac{(l-1)(\mu-\beta)\beta}{\mu}))$ and

$\sum\limits_{i=1}^{l}\ln(1+\frac{i(\mu-\beta)}{k})dx=\int_{1}^{l}\ln(1+\frac{x(\mu-\beta)}{k})dx+O(\ln(1+\frac{l(\mu-\beta)}{k}))$

We integrate the series $\ln(1+x)=x-\frac{1}{2}x^{2}+\frac{1}{3}x^{3}-\cdots$, and take $l=O(k^{2/3})$

[TABLE]

and

[TABLE]

Hence

[TABLE]

or, setting $\sigma^{2}=\frac{k\mu}{(\mu-\beta)^{2}}$

[TABLE]

From Eq.2.13

[TABLE]

and, using Stirling’s formula and the fact that $n_{0}\sim\frac{k}{\mu-\beta}$

[TABLE]

where $\sigma=\sqrt{\frac{\mu k}{(\mu-\beta)^{2}}}$. Thus

[TABLE]

and so

[TABLE]

∎

2.3 Global Limit Theorems

A functional Law of Large Numbers follows directly from Theorem 3.1 in Kurtz (1970 [12]). Likewise, a functional Central Limit Theorem follows from Theorems 3.1 and 3.5 in Kurtz (1971 [13]). We state these theorems here, therefore, without proof.

Write the population size as $n_{k}(t)$, a function of the immigration rate as well as time. Set $n_{k}^{*}=\frac{k}{\mu-\beta}$, the limit of the first moment as $t\to\infty$. Define a new stochastic process for the population size divided by the immigration rate, $Z_{k}(t):=\frac{n_{k}(t)}{k}$. Set $z^{*}=\frac{n_{k}^{*}}{k}=\frac{1}{\mu-\beta}$.

We define the transition function, $f_{k}(\frac{n_{k}}{k},j):=\frac{1}{k}p(n_{k},n_{k}+j)$. Thus,

[TABLE]

Note that $f_{k}(z,j)$ does not, in fact, depend on $k$ and we write simply $f(z,j)$.

Theorem 2.2** (Functional LLN).**

Suppose $\lim\limits_{k\to\infty}Z_{k}(0)=z_{0}$. Then, as $k\to\infty$, $Z_{k}(t)\to Z(t)$ uniformly in probability, where $Z(t)$ is a deterministic process, the solution of

[TABLE]

where

[TABLE]

This has the solution

[TABLE]

Next, define $G_{k}(z):=\displaystyle\sum_{j}j^{2}f_{k}(z,j)=(b+\mu)z+1$. This too does not depend on $k$ and we simply write $G(z)$.

Theorem 2.3** (Functional CLT).**

If $\lim\limits_{k\to\infty}\sqrt{k}\left(Z_{k}(0)-z^{*}\right)=\zeta_{0}$, the processes

[TABLE]

converge weakly in the space of cadlag functions on any finite time interval $[0,T]$ to a Gaussian diffusion $\zeta(t)$ with:

initial value $\zeta(0)=\zeta_{0}$, 2. 2)

mean

[TABLE] 3. 3)

variance

[TABLE]

Suppose, moreover, that $F(z_{0})=0$, i.e., $z_{0}=z^{*}$, the equilibrium point. Then, $Z(t)\equiv z_{0}$ and $\zeta(t)$ is an Ornstein-Uhlenbeck process (OUP) with initial value $\zeta_{0}$, infinitesimal drift

[TABLE]

and infinitesimal variance

[TABLE]

Thus, $\zeta(t)$ is normally distributed with mean

[TABLE]

and variance

[TABLE]

3 Spatial Galton-Watson process with immigration and finite Markov chain spatial dynamics

Let $X=\{x,y,\ldots\}$ be a finite set, and define the following parameters.

1. $\beta(x)$ is the rate of duplication at $x\in X$.
2. $\mu(x)$ is the rate of annihilation at $x\in X$.
3. $a(x,y)$ is the rate of transition $x\to y$.
4. $k(x)$ is the rate of immigration into $x\in X$.

We define $\overrightarrow{n}(t)=\{n(t,x),x\in X\}$, the population at moment $t\geq 0$, with $n(t,x)$ the occupation number of site $x\in X$. Letting $\overrightarrow{\lambda}=\{\lambda_{x}\geq 0,x\in X\}$, we write the Laplace transform of the random vector $\overrightarrow{n}(t)\in\mathbb{R}^{N}$, $N=\mathrm{Card}(X)$ as $u(t,\overrightarrow{\lambda})=E\,e^{-(\overrightarrow{\lambda},\overrightarrow{n}(t))}$.

Now we derive the differential equation for $u(t,\overrightarrow{\lambda})$. Denote the $\sigma$-algebra of events before or including $t$ by $\mathcal{F}_{\leq t}$. Setting $\overrightarrow{\varepsilon}(t,dt))=\overrightarrow{n}(t+dt)-\overrightarrow{n}(t)$

[TABLE]

The conditional distribution of $(\overrightarrow{\lambda},\overrightarrow{\varepsilon})$ under $\mathcal{F}_{\leq t}$ is given by the formulas

1. a)

$P\{(\overrightarrow{\lambda},\overrightarrow{\varepsilon}(t,dt))=\lambda_{x}|\mathcal{F}_{\leq t}\}=n(t,x)\beta(x)dt+k(x)dt$

(the birth of a new offspring at site $x$ or the immigration of a new particle into $x\in X$) 2. b)

$P\{(\overrightarrow{\lambda},\overrightarrow{\varepsilon})=\lambda_{y}|\mathcal{F}_{\leq t}\}=n(t,y)\mu(y)dt$

(the death of a particle at $y\in X$) 3. c)

$P\{(\overrightarrow{\lambda},\overrightarrow{\varepsilon})=\lambda_{x}-\lambda_{z}|\mathcal{F}_{\leq t}\}=n(t,x)a(x,z)dt;\ x,z\in X,\ x\neq z$

(transition of a single particle from $x$ to $z$. Then, $n(t+dt,x)=n(t,x)-1$, $n(t+dt,z)=n(t,z)+1$.) 4. d)

$P\{(\overrightarrow{\lambda},\overrightarrow{\varepsilon})=0|\mathcal{F}_{\leq t}\}=1-(\displaystyle\sum_{x\in X}n(t,x)\beta(x))dt-(\displaystyle\sum_{x\in X}k(x))dt-(\displaystyle\sum_{y\in X}n(t,y)\mu(y))dt-(\displaystyle\sum_{x\neq z}n(t,x)a(x,z))dt$

After substitution of these expressions into 3.1 and elementary transformations we obtain

[TABLE]

But

[TABLE]

I.e., finally

[TABLE]

The initial condition is

[TABLE]

(say, $u(0,\overrightarrow{\lambda})=e^{-(\overrightarrow{\lambda},\mathbf{1})}=e^{\sum_{x\in X}\lambda_{x}}$ for $n(0,x)=1$).

Differentiation of 3.2 and the substitution of $\overrightarrow{\lambda}=0$ leads to the equations for the correlation functions (moments) of the field $n(t,x)$, $x\in X$. Put

[TABLE]

Then

[TABLE]

If $a(x,z)=a(z,x)$ then finally

[TABLE]

Here, $A$ is the generator of a Markov chain $A=[a(x,y)]=A^{*}$.

By differentiating equation 3.2 over the variables $\lambda_{x}$, $x\in X$, one can get the equations for the correlation functions

[TABLE]

where $x_{1},\ldots,x_{m}$ are different points of $X$ and $l_{1},\dots,l_{m}\geq 1$ are integers. Of course $k_{l_{1}\ldots l_{m}}(t,x_{1},\ldots,x_{m})=(-1)^{l_{1}+\dots+l_{m}}\frac{\partial^{l_{1}+\dots+l_{m}}n(t,\overrightarrow{x})}{\partial^{l_{1}}\lambda_{x_{1}}\ldots\partial^{l_{m}}\lambda_{x_{m}}}|_{\overrightarrow{\lambda}=0}$. The corresponding equations will be linear. The central point here is that the factors $(e^{\lambda_{x}-\lambda_{z}}-1)$, $(e^{\lambda_{y}}-1)$, and $(e^{-\lambda_{x}}-1)$ are equal to 0 for $\overrightarrow{\lambda}=0$. As a result, the higher order ($n>l_{1}+\ldots+l_{m}$) correlation functions cannot appear in the equations for {$k_{l_{1}\ldots l_{m}}(\cdot)$, $l_{1}+\ldots+l_{m}=n$}.

Consider, for instance, the correlation function (in fact, matrix valued function)

[TABLE]

The method based on generating functions is typical for the theory of branching processes. In the case of processes with immigration, another, Markovian approach gives new results. Let us start from the simplest case, when there is but one site, i.e., $X=\{x\}$. Then, the process $n(t)$, $t\geq 0$ is a random walk with reflection on the half axis $n\geq 0$.

For a general random walk $y(t)$ on the half axis with reflection in continuous time, we have the following facts. Let the process be given by the generator $G=(g(w,z))$, $w,z\geq 0$, where $a_{w}=g(w,w+1)$, $w\geq 0$; $b_{w}=g(w,w-1)$, $w>0$; $g(w,w)=-(a_{w}+b_{w})$, $w>0$; and $g(0,0)=-a_{0}$ (see Fig. 1).

The random walk is recurrent iff the series

[TABLE]

diverges. It is ergodic (positively recurrent) iff the series

[TABLE]

converges. In the ergodic case, the invariant distribution of the random walk $y(t)$ is given by the formula

[TABLE]

(see [1]).

For our random walk, $n(t)$

[TABLE]

and, for $n\geq 1$

[TABLE]

Proposition 3.1**.**

If $\beta>\mu$ the process $n(t)$ is transient and the population $n(t)$ grows exponentially. 2. 2.

If $\beta=\mu$, $k>0$ the process is not ergodic but rather it is zero-recurrent for $\frac{k}{\beta}\leq 1$ and transient for $\frac{k}{\beta}>1$. 3. 3.

If $\beta<\mu$ the process $n(t)$ is ergodic. The invariant distribution for $\beta<\mu$ is given by

[TABLE]

where $\tilde{S}=\displaystyle\sum_{j=1}^{\infty}\frac{k(k+\beta)\cdots(k+\beta(j-1))}{\mu\cdot 2\mu\cdots j\mu}$.

Proof.

1 and 3 follow from Eqs. 3.3, 3.4, and 3.5. If $\beta=\mu$ (but $k>0$), i.e., in the critical case, the process cannot be ergodic because, setting $\alpha=\frac{k}{\beta}-1$, then $\alpha>-1$ and as $n\to\infty$ $\tilde{S}\sim\displaystyle\sum_{n}n^{\alpha}=+\infty$. The process is zero-recurrent, however, for $0<\frac{k}{\beta}\leq 1$. In fact, for $\beta=\mu$

[TABLE]

and the series in Eq. 3.4 diverges if $0<\frac{k}{\beta}\leq 1$. If, however, $k>\beta$ the the series converges and the process $n(t)$ is transient. ∎

Consider, now, the general case of the finite space $X$. Let $N=\mathrm{Card}\,X$ and $\overrightarrow{n}(t)$ be the vector of the occupation numbers. The process $\overrightarrow{n}(t)$, $t\geq 0$ is the random walk on $(\mathbb{Z}_{+}^{1})^{N}=\{0,1,...)^{N}$ with continuous time. The generator of this random walk was already described when we calculated the Laplace transform $u(t,\overrightarrow{\lambda})=E\,e^{-(\overrightarrow{\lambda},\overrightarrow{n}(t))}$. If at the moment $t$ we have the configuration $\overrightarrow{n}(t)=\{n(t,x),x\in X\}$, then, for the interval $(t,t+dt)$ only the following events (up to terms of order$(dt)^{2}$) can happen:

1. a)

the birth of a new particle at the site $x_{0}\in X$, with corresponding probability $n(t,x_{0})\beta(x_{0})dt+k(x_{0})dt$. In this case we have the transition

[TABLE] 2. b)

the death of one particle at the site $x_{0}\in X$. This has corresponding probability $\mu(x_{0})n(t,x_{0})dt$ and the transition

[TABLE]

(Of course, here $n(t,x_{0})\geq 1$, otherwise $\mu(x_{0})n(t,x_{0})dt=0$). 3. c)

the transfer of one particle from site $x_{0}$ to site $y_{0}\in X$ (jump from $x_{0}$ to $y_{0}$), i.e., the transition

[TABLE]

with probability $n(t,x_{0})a(x_{0},y_{0})dt$ for $n(t,x_{0})\geq 1$.

The following theorem gives sufficient conditions for the ergodicity of the process $\overrightarrow{n}(t)$.

Theorem 3.2**.**

Assume that for some constants $\delta>0$, $A>0$ and any $x\in X$

[TABLE]

Then, the process $\overrightarrow{n}(t)$ is an ergodic Markov chain and the invariant measure of this process has exponential moments, i.e., $E\,e^{(\overrightarrow{\lambda},\overrightarrow{n}(t))}\leq c_{0}<\infty$ if $|\overrightarrow{\lambda}|\leq\lambda_{0}$ for appropriate (small) $\lambda_{0}>0$.

Proof.

We take on $(\mathbb{Z}_{+}^{1})^{N}=\{0,1,...)^{N}$ as a Lyapunov function

[TABLE]

Then, with $G$ the generator of the process, $Gf(\overrightarrow{n})\leq 0$ for large enough $(\overrightarrow{n},\overrightarrow{\mathbf{1}})=\|\overrightarrow{n}\|_{1}$. In fact

[TABLE]

(The terms concerning transitions of the particles between sites make no contribution: $1-1=0$.) ∎

If $\beta(x)\equiv\beta<\mu\equiv\mu(x)$ and $k(x)\equiv k$ then $(\overrightarrow{n},\overrightarrow{\mathbf{1}})$, i.e., the total number of the particles in the phase space $X$ is also a Galton Watson process with immigration and the rates of transition shown in Fig. 2.

If $t\to\infty$ this process has a limiting distribution with invariant measure (in which $Nk$ replaces $k$). That is

[TABLE]

4 Branching process with migration and immigration

We now consider our process with birth, death, migration, and immigration on a countable space, specifically the lattice $\mathbb{Z}^{d}$. As in the other models, we have $\beta>0$, the rate of duplication at $x\in Z^{d}$; $\mu>0$, the rate of death; and $k>0$, the rate of immigration. Here, we add migration of the particles with rate $\kappa>0$ and probability kernel $a(z)$, $z\in\mathbb{Z}^{d}$, $z\neq 0$, $a(z)=a(-z)$, $\sum\limits_{z\neq 0}a(z)=1$. That is, a particle jumps from site $x$ to $x+z$ with probability $\kappa a(z)dt$. Here we put $\kappa=1$ to simplify the notation.

For $n(t,x)$ the number of particles at $x$ at time $t$, the forward equation for this process is given by $n(t+dt,x)=n(t,x)+\xi(dt,x)$, where

[TABLE]

Note that $\xi(dt,x)$ is independent on $\mathcal{F}_{\leqslant t}$ (the $\sigma$-algebra of events before or including $t$) and

1. a)

$E[\xi(dt,x)|\mathcal{F}_{\leqslant t}]=n(t,x)(\beta-\mu-1)dt+kdt+\sum\limits_{z\neq 0}a(z)n(t,x+z)dt$. 2. b)

$E[\xi^{2}(dt,x)|\mathcal{F}_{\leqslant t}]=n(t,x)(\beta+\mu+1)dt+kdt+\sum\limits_{z\neq 0}a(z)n(t,x+z)dt$. 3. c)

$E[\xi(dt,x)\xi(dt,y)|\mathcal{F}_{\leqslant t}]=a(x-y)n(t,x)dt+a(y-x)n(t,y)dt$.

A single particle jumps from $x$ to $y$ or from $y$ to $x$. Other possibilities have probability $O((dt)^{2})\approx 0$. Here, of course, $x\neq y$. 4. d)

If $x\neq y$, $y\neq z$, and $x\neq z$, then $E[\xi(dt,x)\xi(dt,y)\xi(dt,z)]=0$.

We will not use property d) in this paper but it is crucial for the analysis of moments of order greater or equal to $3$.

From here on, we concentrate on the first two moments.

4.1 First moment

Due to the fact that $\beta<\mu$, the system has a short memory, and we can calculate all the moments under the condition that $n(0,x)$, $x\in\mathbb{Z}^{d}$, is a system of independent and identically distributed random variables with expectation $\frac{k}{\mu-\beta}$. We will select Poissonian random variables with parameter $\lambda=\frac{k}{\mu-\beta}$. Then, $m_{1}(t,x)=\frac{k}{\mu-\beta}$, $t\geqslant 0$, $x\in\mathbb{Z}^{d}$, and, as a result, $\mathcal{L}_{a}m_{1}(t,x)=0$. Setting $m_{1}(t,x)=E[n(t,x)]$, we have

[TABLE]

Defining the operator $\mathcal{L}_{a}(f(t,x))=\sum\limits_{z\neq 0}a(z)[f(t,x+z)-f(t,x)]$, then, from Eq. 4.2 we get the differential equation

[TABLE]

Because of spatial homogeneity, $\mathcal{L}_{a}m_{1}(t,x)=0$, giving

[TABLE]

which has the solution

[TABLE]

Thus, if $\beta\geq\mu$, $m_{1}(t,x)\rightarrow\infty$, and if $\mu>\beta$,

[TABLE]

4.2 Second moment

We derive differential equations for the second correlation function $m_{2}(t,x,y)$ for $x=y$ and $x\neq y$ separately, then combine them and use a Fourier transform to prove a useful result concerning the covariance.

1. I.

$x=y$

[TABLE]

Denote $\mathcal{L}_{ax}m_{2}(t,x,y)=\sum\limits_{z\neq 0}a(z)(m_{2}(t,x+z,y)-m_{2}(t,x,y))$.

From this follows the differential equation

[TABLE] 2. II.

$x\neq y$

Because only one event can happen during $dt$

[TABLE]

while the probability that one particle jumps from $y$ to $x$ is

[TABLE]

and the probability that one particle jumps from $x$ to $y$ is

[TABLE]

Then, similar to above

[TABLE]

The resulting differential equation is

[TABLE]

That is

[TABLE]

Because, for fixed $t$, $n(t,x)$ is homogeneous in space, we can write $m_{2}(t,x,y)=m_{2}(t,x-y)=m_{2}(t,u)$. Then, we can condense the two cases into a single differential equation

[TABLE]

Here $u=x-y\neq 0$ and $a(0)=0$.

We can partition $m_{2}(t,u)$ into $m_{2}(t,u)=m_{21}+m_{22}$, where the solution for $m_{21}$ depends on time but not position and the solution for $m_{22}$ depends on position but not time. Thus, $\mathcal{L}_{au}m_{21}=0$ and $m_{21}$ corresponds to the source $\frac{2k^{2}}{\mu-\beta}$, which gives

[TABLE]

As $t\rightarrow\infty$, $m_{21}\rightarrow\bar{M_{2}}=m_{1}^{2}(t,x)=\frac{k^{2}}{(\mu-\beta)^{2}}$.

For the second part, $m_{22}$ , $\frac{\partial m_{22}}{\partial t}=0$, i.e.

[TABLE]

As $t\rightarrow\infty$, $m_{22}\rightarrow\tilde{M}_{2}$. $\tilde{M}_{2}$ is the limiting correlation function for the particle field $n(t,x)$, $t\rightarrow\infty$. It is the solution of the “elliptic” problem

[TABLE]

Applying the Fourier transform $\widehat{\tilde{M_{2}}}(\theta)=\sum\limits_{u\in Z^{d}}\tilde{M_{2}}(u)e^{i(\theta,u)}$, $\theta\in T^{d}=[-\pi,\pi]^{d}$,

we obtain

[TABLE]

We have proved the following result.

Theorem 4.1**.**

*If $t\rightarrow\infty$, then $\mathrm{Cov}(n(t,x),n(t,y))=E[n(t,x)n(t,y)]-E[n(t,x)]E[n(t,y)]\\ =\-m_{2}(t,x,y)-m_{1}(t,x)m_{1}(t,y)$, tends to $\tilde{M_{2}}(x-y)=\tilde{M_{2}}(u)\in L^{2}(Z^{d})$

The Fourier transform of $\tilde{M_{2}}(\cdot)$ is equal to*

[TABLE]

where $c_{1}=\frac{k}{\mu-\beta}$, $c_{2}=\frac{k}{\mu-\beta}$, $c_{3}=\mu-\beta$

Let’s compare our results with the corresponding results for the critical contact model [2] (where $k=0$, $\mu=\beta$). In the last case, the limiting distribution for the field $n(t,x)$, $t\geqslant 0$, $x\in\mathbb{Z}^{d}$, exists if and only if the underlying random walk with generator $\mathcal{L}_{a}$ is transient. In the recurrent case, we have the phenomenon of clusterization. The limiting correlation function is always slowly decreasing (like the Green kernel of $\mathcal{L}_{a}$).

In the presence of immigration, the situation is much better: the limiting correlation function always exists and we believe that the same is true for all moments. The decay of $\tilde{M}_{2}(u)$ depends on the smoothness of $\hat{a}(\theta)$. Under minimal regularity conditions, correlations have the same order of decay as $a(z)$, $z\rightarrow\infty$. For instance, if $a(z)$ is finitely supported or exponentially decreasing, the correlation also has an exponential decay. If $a(z)$ has power decay, then the same is true for correlation $\tilde{M}_{2}(u)$, $u\rightarrow\infty$.

5 Processes in a Random Environment

The final four models involve a random environment. Two are Galton-Watson models with immigration and lack a spatial component. In the first, the parameters are random functions of the population size; in the second, they are random functions of a Markov chain on a finite space. The last two models are spatial and feature immigration, migration, and, most importantly, a random environment in space, still stationary in time for the third but not stationary in time for the fourth.

5.1 Galton-Watson processes with immigration in random environments

5.1.1 Galton-Watson process with immigration in random environment based on population size

Assume that rates of mortality $\mu(\cdot)$, duplication $\beta(\cdot)$, and immigration $k(\cdot)$ are random functions of the volume of the population $x\geq 0$. Namely, the random vectors $(\mu,\beta,k)(x,\omega)$ are i.i.d on the underlying probability space $(\Omega_{e},\mathcal{F}_{e},P_{e})$ ($e$: environment).

The Galton-Watson Process is ergodic ($P_{e}$-a.s) if and only if the random series

[TABLE]

Theorem 5.1**.**

Assume that the random variables $\beta(x,\omega)$, $\mu(x,\omega)$, $k(x,\omega)$ are bounded from above and below by the positive constants $C^{\pm}$: $0<C^{-}\leq\beta(x,\omega)\leq C^{+}<\infty$. Then, the process $n(t,\omega_{e})$ is ergodic $P_{e}$-a.s. if and only if $\langle\ln\frac{\beta(x,\omega)}{\mu(x,\omega)}\rangle=\langle\ln\beta(\cdot)\rangle-\langle\ln(\mu(\cdot))\rangle<0$

Proof.

It is sufficient to note that

[TABLE]

∎

It follows from the strong LLN that the series diverges exponentially fast for $\langle\ln\beta(\cdot)\rangle-\langle\ln\mu(\cdot)\rangle>0$; it converges like a decreasing geometric progression for $\langle\ln\beta(\cdot)\rangle-\langle\ln\mu(\cdot)\rangle<0$; and it is divergent if $\langle\ln\beta(\cdot)\rangle=\langle\ln\mu(\cdot)\rangle$. It diverges even when $\beta(x,\omega_{e})=\mu(x,\omega_{e})$ due to the presence of $k^{-}\geq C^{-}>0.$

Note that $ES<\infty$ if and only if $\langle\frac{\lambda(x-1)}{\mu(x)}\rangle=\langle\lambda\rangle\langle\frac{1}{\mu}\rangle<\infty$, i.e., the fluctuations of $S$, even in the case of convergence, can be very high.

5.1.2 Random non-stationary(time dependent) environment

Assume that $k(t)$ and $\Delta=(\mu-\beta)(t)$ are stationary random processes on $(\Omega_{m},P_{m})$ and that $k(t)$ is independent of $\Delta$. For a fixed environment, i.e., fixed $k(\cdot)$ and $\Delta(\cdot)$, the equation for the first moment takes the form

[TABLE]

Then

[TABLE]

Assume that $\frac{1}{\delta}\geqslant\Delta(\cdot)\geqslant\delta>0$, $\frac{1}{\delta}\geqslant k(\cdot)\geqslant\delta>0$. Then

[TABLE]

Thus, for large $t$, the process $m_{1}(t,\omega_{m})$ is exponentially close to the stationary process

[TABLE]

Assume now that $k(t)$ and $\Delta(s)$ are independent stationary processes and $-\Delta(t)=V(x(t))$, where $x(t)$, $t\geqslant 0$, is a Markov Chain with continuous time and symmetric geometry on the finite set $X$. (One can also consider $x(t)$, $t\geqslant 0$, as a diffusion process on a compact Riemannian manifold with Laplace-Beltrami generator $\Delta$.) Let

[TABLE]

Then

[TABLE]

The operator $\mathcal{L}$ is symmetric in $L^{2}(x)$ with dot product $(f,g)=\sum\limits_{x\in X}f(x)\bar{g(x)}$. Thus, $H=\mathcal{L}+V$ is also symmetric and has real spectrum $0>-\delta\geqslant\lambda_{0}>\lambda_{1}\geqslant\cdots$ with orthonormal eigenfunctions $\psi_{0}(x)>0$,$\psi_{1}(x)>0$, $\cdots$ Iinequality $\lambda_{0}\leqslant\delta<0$ follows from our assumption on $\Delta(\cdot)$.

The solution of equation 5.1 is given by

[TABLE]

Now, we can calculate $<\tilde{m_{1}}(t,x,\omega_{m})>$.

[TABLE]

Here, $\pi(x)=\frac{1}{N}=\frac{\mathbbm{1}(x)}{N}$ is the invariant distribution of $x_{s}$. Then

[TABLE]

5.1.3 Galton-Watson process with immigration in random environment given by Markov chain

Let $x(t)$ be an ergodic MCh on the finite space $X$ and let $\beta(x),\mu(x),k(x)$, the rates of duplication, annihilation, and immigration, be functions from $X$ to $R^{+}$, and, therefore, functions of $t$ and $\omega_{e}$. The process $(n(t),x(t))$ is a Markov chain on $\mathbb{Z}_{+}^{1}\times X$.

Let $a(x,y)$, $x,y\in X$, $a(x,y)\geq 0$, $\sum\limits_{y\in X}a(x,y)=1$ for all $x\in X$, be the transition function for $x(t)$. Consider $E_{(n,x)}f(n(t),x(t))=u(t,(n,x))$. Then

[TABLE]

We obtain the backward Kolmogorov equation

[TABLE]

Example. Two-state random environment.

Here, $x(t)$ indicates which one of two possible states, $\{1,2\}$ the process is in at time $t$. The birth, mortality, and immigration rates are different for each state: $\beta_{1}$ and $\beta_{2}$, $\mu_{1}$ and $\mu_{2}$, and $k_{1}$ and $k_{2}$. For a process in state 1, at any time the rate of switching to state 2 is $\alpha_{1}$, with $\alpha_{2}$ the rate of the reverse switch. This creates the two-state random environment. Let $G$ be the generator for the process, as diagrammed in Figure 3.

The following theorem gives sufficient conditions for the ergodicity of the process $(n(t),x(t))$.

Theorem 5.2**.**

Assume that for some constants $\delta>0$ and $A>0$

[TABLE]

Then, the process $(n(t),x(t))$ is an ergodic Markov chain and the invariant measure of this process has exponential moments, i.e., $E\,e^{\lambda n(t)}\leq c_{0}<\infty$ if $\lambda\leq\lambda_{0}$ for appropriate (small) $\lambda_{0}>0$.

Proof.

We take as a Lyapunov function $f(n,x)=n$.

Then, $Gf(n(t),x(t))=(\beta_{x}-\mu_{x})n+k_{x}$. So for sufficiently large $n$, specifically $n>\frac{A}{\delta}$, we have $Gf\leq 0$. ∎

5.2 Models with immigration and migration in a random environment

For this most general case, we have migration and a non-stationary environment in space and time. The rates of duplication, mortality, and immigration at time $t$ and position $x\in\mathbb{Z}^{d}$ are given by $\beta(t,x)$, $\mu(t,x)$, and $k(t,x)$. As in the above models, immigration is uninfluenced by the presence of other particles; also set $\delta_{1}\leq k(t,x)\leq\delta_{2}$, $0<\delta_{1}<\delta_{2}<\infty$. The rate of migration is given by $\kappa$, with the process governed by the probability kernel $a(z)$, the rate of transition from $x$ to $x+z$, $z\in Z^{d}$.

If $n(t,x)$ is the number of particles at $x\in Z^{d}$ at time $t$, $n(t+dt,x)=n(t,x)+\xi(t,x)$, where

[TABLE]

For the first moment, $m_{1}(t,x)=E[n(t,x)]$, we can write

[TABLE]

and so, defining, as above, $\mathcal{L}_{a}(f(t,x))=\sum\limits_{z\neq 0}a(z)[f(t,x+z)-f(t,x)]$, we obtain

[TABLE]

We consider two cases. The first is where the duplication and mortality rates are equal, $\beta(t,x)=\mu(t,x)$. Because of the immigration rate bounded above 0, we find that the expected population size at each site tends to infinity. In the second case, to simplify, we consider $\beta(t,x)$ and $\mu(t,x)$ to be stationary in time, and assume the mortality rate to be greater than the duplication rate everywhere by at least a minimal amount. Here, we show that the interplay between the excess mortality and the positive immigration results in a finite positive expected population size at each site.

5.2.1 Case I

If $\beta(t,x)=\mu(t,x)$

[TABLE]

Taking Fourier transforms,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Taking the inverse Fourier transform,

[TABLE]

where

[TABLE]

As $t\rightarrow\infty$, $\delta_{1}t\rightarrow\infty$. Thus, when the birth rate equals the death rate, the expected population at each site $x\in\mathbb{Z}^{d}$ will go to infinity as $t\rightarrow\infty$.

5.2.2 Case II

Here, $\beta(t,x)\neq\mu(t,x)$. For simplification we assume that only immigration, $k(t,x)$, is not stationary in time. In other words, we us assume that the duplication and mortality rates are stationary in time and depend only on position: $\beta(t,x)=\beta(x)$, $\mu(t,x)=\mu(x)$ and $\mu(x)-\beta(x)\geqslant\delta_{1}>0$. From Eq. 5.3, we get

[TABLE]

This has the solution

[TABLE]

where $q(t-s,x,y)$ is the solution for

[TABLE]

By the Feynman-Kac formula,

[TABLE]

where

[TABLE]

Finally

[TABLE]

Thus, when $\mu(x)-\beta(x)>0$, $\lim\limits_{t\rightarrow\infty}m_{1}(t,x)$ is bounded by [math] and $\frac{\lVert k\rVert_{\infty}}{\delta_{1}}$, so this limit exists and is finite.