A Population with Immigration and Density-Dependent Deaths
In a population model without immigration, the state zero is often absorbing: once the population becomes extinct, it cannot return. Immigration changes that mechanism. Even when the population is empty, an arrival can restart it. To describe the resulting long-run population sizes, we will derive stationary probabilities from the balance of transitions between each pair of adjacent population sizes.
Work on a probability space \((\Omega,\mathcal{F},\mathbb{P})\), and let \((X_t)_{t\ge0}\) be a continuous-time Markov chain with state space \(S=\mathbb{N}_0=\{0,1,2,\ldots\}\), natural filtration \((\mathcal{F}_t)_{t\ge0}\), and rates specified below. A transition from \(n\) to \(n+1\) is a birth or an immigrant; a transition from \(n\) to \(n-1\) is a death. The total rate of increasing the population from \(n\) is \(\lambda_n=a+bn\), where \(a>0\) is the immigration rate and \(b\ge0\) is the per-capita birth rate. The death rate is \(\mu_n=n[d+c(n-1)]\) for \(n\ge1\), with \(\mu_0=0\), \(d>0\), and \(c>0\). The term involving \(c\) makes deaths increase more rapidly than linearly as the population grows.
The adjacent-state balance equation equates the stationary flow from \(n\) to \(n+1\) with the flow back from \(n+1\) to \(n\): \(\pi_n\lambda_n=\pi_{n+1}\mu_{n+1}\). This equation suggests all the stationary probabilities once \(\pi_0\) is chosen. The remaining question is whether the resulting nonnegative weights have a finite sum, so that they can be normalized to a probability distribution.
Proof. First, the weights satisfy \(w_n\lambda_n=w_{n+1}\mu_{n+1}\) for every \(n\ge0\), directly from the product defining \(w_{n+1}\). They are positive because \(a,d,c>0\). To check that they can be normalized, their successive ratio is $$ \frac{w_n}{w_{n-1}}=\frac{a+b(n-1)}{n[d+c(n-1)]},\qquad n\ge1. $$ As \(n\) grows, the numerator is at most a constant times \(n\), while the denominator grows at least as a positive constant times \(n^2\). Consequently the ratio tends to zero. In particular, it is at most \(1/2\) for all sufficiently large \(n\), so the tail of the series of weights is bounded by a geometric series. Thus \(Z<\infty\).
For completeness, the rates define a non-explosive chain, meaning the chain makes only finitely many jumps in any bounded time interval. Births can be coupled below a pure-birth population with rates \(a+bn\): the latter has every possible upward jump of the original chain, while deaths in the original chain only reduce its population. The pure-birth process is non-explosive; for example, its expected size stopped on first reaching level \(M\) is bounded on each fixed time interval by the solution of the linear differential inequality \(u'(t)\le a+bu(t)\), independently of \(M\). Markov's inequality then makes the probability of reaching \(M\) by that time tend to zero as \(M\) tends to infinity. This also rules out explosion in the original process.
Now normalize the weights, so \(\pi_n=w_n/Z\), and note that they satisfy \(\pi_nq_{n,n+1}=\pi_{n+1}q_{n+1,n}\). To see why this gives stationarity here, consider any finite-jump path from state \(i\) to state \(j\) over a fixed time interval. Its probability density includes the product of its jump rates and the exponential holding-time factors. Reversing the path changes the jump-rate product by the factor \(\pi_j/\pi_i\), by applying the adjacent-state identity at each jump; the holding-time factors are unchanged because the same states and holding durations occur in reverse order. Integrating over jump times and summing over paths gives \(\pi_iP_{ij}(t)=\pi_jP_{ji}(t)\), where \(P_{ij}(t)=\mathbb{P}_i(X_t=j)\). Non-explosion ensures these paths account for the transition probabilities. Summing over \(i\) yields \(\sum_i\pi_iP_{ij}(t)=\pi_j\sum_iP_{ji}(t)=\pi_j\), since the transition probabilities out of state \(j\) sum to one. This proves stationarity.
Finally, any stationary distribution must obey the same adjacent-state equations. One way to see the necessity is to balance the total stationary flow across the cut \(\{0,\ldots,n\}\): the only transitions crossing its boundary are \(n\) to \(n+1\) and \(n+1\) to \(n\), so their flows must agree. These equations recursively determine every probability from \(\pi_0\). Normalization then forces \(\pi_0=1/Z\), proving uniqueness. \(\square\)
Calculating and Interpreting the Stationary Probabilities
The product formula is the main calculation: start with \(w_0=1\), multiply by the ratio \(\lambda_{n-1}/\mu_n\) to obtain each next weight, and divide all weights by their sum. In general, that sum need not simplify to an elementary expression. The examples below use rates for which it can be evaluated exactly.
Worked Example: Immigration Rate One
Take \(a=1\), \(b=1\), \(d=2\), and \(c=1\). Then \(\lambda_n=1+n\) and \(\mu_n=n(n+1)\) for \(n\ge1\). The generator row at zero sums to \(-1+1=0\); for \(n\ge1\), its three possibly nonzero entries sum to \((n+1)+n(n+1)-[(n+1)+n(n+1)]=0\). Thus the rates give a valid generator.
For \(k\ge1\), the ratio in the weight formula is $$ \frac{a+b(k-1)}{k[d+c(k-1)]} =\frac{1+(k-1)}{k[2+(k-1)]} =\frac{k}{k(k+1)} =\frac{1}{k+1}. $$ Therefore \(w_n=\prod_{k=1}^{n}1/(k+1)=1/(n+1)!\), including \(w_0=1\). The normalizing sum is \(Z=\sum_{n=0}^{\infty}1/(n+1)!=e-1\), so \(\pi_n=1/[(e-1)(n+1)!]\).
Check: For any \(n\ge0\), \(\pi_n\lambda_n=(n+1)/[(e-1)(n+1)!]=1/[(e-1)n!]\). Also, \(\pi_{n+1}\mu_{n+1}=(n+1)(n+2)/[(e-1)(n+2)!]=1/[(e-1)n!]\). The adjacent-state flows agree.
Worked Example: A Higher Immigration Rate
Keep \(b=1\), \(d=2\), and \(c=1\), but set \(a=2\). Here \(\lambda_n=2+n\) and \(\mu_n=n(n+1)\). The generator row at zero sums to \(-2+2=0\); for \(n\ge1\), the entries sum to \((n+2)+n(n+1)-[(n+2)+n(n+1)]=0\).
For \(k\ge1\), the weight ratio is $$ \frac{2+(k-1)}{k[2+(k-1)]} =\frac{k+1}{k(k+1)} =\frac{1}{k}. $$ Hence \(w_n=\prod_{k=1}^n1/k=1/n!\), and \(Z=\sum_{n=0}^{\infty}1/n!=e\). The stationary distribution is \(\pi_n=1/(e\,n!)\).
Check: The adjacent-state flows are \(\pi_n\lambda_n=(n+2)/(e\,n!)\) and \(\pi_{n+1}\mu_{n+1}=(n+1)(n+2)/(e\,(n+1)!)=(n+2)/(e\,n!)\). They agree, and \(\sum_{n=0}^{\infty}\pi_n=e/e=1\).
In these two examples, the stationary probability of an empty population falls from \(1/(e-1)\) to \(1/e\) when the immigration rate increases from one to two. This is not just a feature of these selected rates. For fixed \(b,d,c\), each weight \(w_n\) for \(n\ge1\) increases strictly with \(a\): the factor for \(k=1\) is \(a/d\), and all remaining factors are nondecreasing in \(a\). Therefore \(Z=\sum_n w_n\) increases strictly with \(a\), and \(\pi_0=1/Z\) decreases strictly. This comparison concerns the stationary probability at one state; it does not by itself assert a particular ordering for every other population probability.
Worked Example: Removing Immigration
Set \(a=0\), \(b=1\), \(d=2\), and \(c=1\). Now \(\lambda_n=n\), so \(\lambda_0=0\) and state zero is absorbing. For \(n\ge1\), \(\mu_n=n(n+1)\). The generator row at zero sums to \(0\), while each row \(n\ge1\) sums to \(n+n(n+1)-[n+n(n+1)]=0\).
The distribution concentrated at zero, \(\pi_0=1\) and \(\pi_n=0\) for \(n\ge1\), is stationary because no transition leaves zero. It is the only stationary probability distribution: balance across the cut \(\{0\}\) requires \(\pi_0\lambda_0=\pi_1\mu_1\), so \(0=2\pi_1\) and \(\pi_1=0\). Balance across \(\{0,1\}\) then gives \(0=\pi_2\mu_2\), hence \(\pi_2=0\); repeating this argument forces every \(\pi_n\) for \(n\ge1\) to be zero. The total-mass condition gives \(\pi_0=1\).
Check: The only positive mass is at an absorbing state, so evolving this distribution leaves it unchanged. By contrast, when \(a>0\), zero has a positive rate \(a\) of transition to one and cannot be absorbing.
What Stationary Persistence Means
With positive immigration, every state communicates with every other: consecutive upward transitions have positive rates because \(\lambda_n\ge a>0\), and consecutive downward transitions have positive rates because \(\mu_n>0\) for \(n\ge1\). The stationary distribution assigns positive probability to every population size, including zero and positive sizes. Thus immigration makes the empty state non-absorbing and gives the population a stationary distribution with positive mass away from zero.
This is a precise but limited sense of persistence. The population may still visit zero; immigration gives it a way to leave zero again. A stationary probability is an invariant distribution, not by definition a statement that the distribution at time \(t\) converges to it from every starting state. No such limiting claim is needed for the balance calculation here. It is also important not to infer survival forever from a positive stationary mean or from \(\pi_0<1\).
Common mistake. A tempting but incorrect calculation is to set \(\pi_n\lambda_n=\pi_n\mu_n\), comparing the upward and downward rates at the same state. That is not the flow across the boundary between \(n\) and \(n+1\): the downward flow across that boundary starts at \(n+1\). The correct equation is \(\pi_n\lambda_n=\pi_{n+1}\mu_{n+1}\). For the first example, at the boundary between zero and one, the correct equation is \(\pi_0\cdot1=\pi_1\cdot2\), giving \(\pi_1=\pi_0/2\); the incorrect same-state equation would instead demand \(\pi_0\cdot1=\pi_0\cdot0\), which is impossible for a probability with \(\pi_0>0\).
Check Your Understanding
Use the rate definitions and balance equations to answer each question. In each case, include the relevant boundary condition.
- For the model with \(a>0\), what are the upward and downward rates across the boundary between states \(n\) and \(n+1\)?
- With \(a=1,b=1,d=2,c=1\), compute \(w_1,w_2\), and verify the adjacent-state balance equation at \(n=1\).
- Why does the density-dependent term ensure that the stationary weights have a finite sum for the stated parameter conditions?
- For fixed \(b,d,c\), explain why increasing \(a\) lowers the stationary probability of zero.
- When \(a=0\), explain why the stationary distribution is concentrated at zero and why this does not describe a population that can recover after reaching zero.