# On Dirichlet eigenvectors for neutral two-dimensional Markov chains

Nicolas Champagnat, Persi Diaconis, Laurent Miclo

► To cite this version:

Nicolas Champagnat, Persi Diaconis, Laurent Miclo. On Dirichlet eigenvectors for neutral two-dimensional Markov chains. Electronic Journal of Probability, 2012, 17 (63), pp.1-41. 10.1214/EJP.v17-1830. hal-00672938

HAL Id: hal-00672938

https://inria.hal.science/hal-00672938v1

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On Dirichlet eigenvectors for neutral two-dimensional Markov chains

Nicolas Champagnat,\* Persi Diaconis,† Laurent Miclo‡

## Abstract

We consider a general class of discrete, two-dimensional Markov chains modeling the dynamics of a population with two types, without mutation or immigration, and neutral in the sense that type has no influence on each individual's birth or death parameters. We prove that all the eigenvectors of the corresponding transition matrix or infinitesimal generator $\Pi$ can be expressed as the product of explicit “universal” polynomials of two variables, depending on each type's size but not on the specific transitions of the dynamics, and functions depending only on the total population size. We also prove that all the Dirichlet eigenvectors of $\Pi$ on subdomains of the form $\{(i,j) \in \mathbb{N}^2 : i + j \ge N\}$ for some $N \ge 2$ have the same decomposition. We then deduce that all the corresponding Dirichlet eigenvalues are ordered in a specific way and each of them is related to the greatest eigenvalue associated to eigenvectors admitting one specific “universal” polynomial as factor. As an application, we study the quasistationary behavior of finite, two-dimensional Markov chains such that 0 is an absorbing state for each component of the process. In particular, we prove that coexistence is never possible conditionally on non-extinction in a population close to neutrality.

<footnote>\*TOSCA project-team, INRIA Nancy – Grand Est, IECN, UMR 7502, Nancy-Université, Campus Scientifique, B.P. 70239, 54506 Vandœuvre-lès-Nancy Cedex, France, E-mail: Nicolas.Champagnat@inria.fr</footnote>

<footnote>†Department of Mathematic and Statistics, Sequoia Hall, 390 Serra Mall, Stanford University, Stanford, CA 94305-4065s</footnote>

<footnote>‡Institut de Mathématiques de Toulouse, Université Paul Sabatier, 118 route de Narbonne, F-31062 Toulouse Cedex 9, France</footnote>
MSC 2000 subject classification: Primary: 60J10, 60J27; secondary: 15A18, 39A14, 47N30, 92D25.

Keywords and phrases: Hahn polynomials; two-dimensional difference equation; neutral Markov chain; multitype population dynamics; Dirichlet eigenvector; Dirichlet eigenvalue; quasi-stationary distribution; Yaglom limit; coexistence.

# 1 Introduction

This paper studies spectral properties of two-dimensional discrete Markov processes in continuous and discrete time, having the neutrality property, in the sense of population genetics (see e.g. [19]). Considering two populations in interaction, corresponding to two different types of individuals (typically a mutant type and a resident type), one says that the types are neutral (or the mutant type is neutral w.r.t. the resident type, or more simply the mutation is neutral) if individuals of both types are indistinguishable in terms of the total population dynamics. In other words, the mutant population has no selective advantage (or disadvantage) with respect to the rest of the population.

We consider two-dimensional Markov processes $(X_n, Y_n)_{n \in \mathcal{T}}$ (where $\mathcal{T} = \mathbb{Z}_+ := \{0, 1, \dots\}$ or $\mathbb{R}_+$) with values in $\mathbb{R}_+^2$ or $\mathbb{Z}_+^2$. The notion of neutrality we will consider can be formally defined as follows. In the case of continuous time, we assume that the birth and death rates per individual do not depend on the type of the individual, and depend only on the total population size. In the case of discrete time, we assume that the transition probabilities can be constructed by first determining the number of births or deaths in the next time step, with distribution depending only on the total population size, and second by choosing uniformly at random the individuals concerned by the birth or death events (regardless of their types). As a consequence, the process $(X_n, Y_n)_{n \in \mathcal{T}}$ is such that $Z_n = X_n + Y_n$ is a Markov process. In particular, the law of the process $Z$ depends on $Z_0$, but not on $X_0$ or $Y_0$. Note that other notions of neutrality different from the one considered here can be relevant in biology and can satisfy the last property (cf. e.g. [1]).

If the process $Z$ is a birth and death continuous-time chain, the class of neutral processes we consider is the following: the birth and death rates of the Markov process $(Z_t)_{t \in \mathbb{R}_+}$ when $Z$ is in state $k \ge 0$ are of the form $k\lambda_k$ and $k\mu_k$, respectively. Note that $0$ is an absorbing state for $Z$. In other words, the process $(X_t, Y_t)_{t \in \mathbb{R}_+}$ is the birth and death process where both types of
individuals have birth and death rates per individual $\lambda_k$ and $\mu_k$, when the total population has $k$ individuals. This leads to the following transition rates for the Markov process $(X_t, Y_t)_{t \in \mathbb{R}_+}$: for all $(i, j) \in \mathbb{Z}_+^2$,

<table><tr><td>from (i, j) to (i + 1, j)</td><td>with rate</td><td>i &lambda;<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i, j + 1)</td><td>with rate</td><td>j &lambda;<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i &minus; 1, j)</td><td>with rate</td><td>i &mu;<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i, j &minus; 1)</td><td>with rate</td><td>j &mu;<sub>i+j</sub>.</td></tr></table>

Note that the sets $\{0\} \times \mathbb{Z}_+$, $\mathbb{Z}_+ \times \{0\}$ and $\{(0, 0)\}$ are absorbing for this process. In other words, we only consider neutral two-dimensional processes *without mutation and immigration*.

In the case of discrete time, we consider two-dimensional birth and death processes constructed in a similar way: assume that the birth and death probabilities of the process $(Z_n)_{n \in \mathbb{Z}_+}$ when in state $k$ are $p_k$ and $q_k$, respectively, with $p_k + q_k \le 1$. Then, when a birth or a death occurs in the population, the individual to which this event applies is chosen uniformly at random in the population. This leads to the transition probabilities

<table><tr><td>from (i, j) to (i + 1, j)</td><td>with probability</td><td>i / (i+j) p<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i, j + 1)</td><td>with probability</td><td>j / (i+j) p<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i &minus; 1, j)</td><td>with probability</td><td>i / (i+j) q<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i, j &minus; 1)</td><td>with probability</td><td>j / (i+j) q<sub>i+j</sub></td></tr><tr><td>from (i, j) to (i, j)</td><td>with probability</td><td>r<sub>k</sub></td></tr></table>

where $r_k := 1 - p_k - q_k$. Note that this construction requires assuming that $r_0 = 1$ (i.e. that 0 is absorbing for $Z$).

In [15], Karlin and McGregor studied two families of neutral multitype population processes (branching processes and Moran model), but only in the case of nonzero mutation or immigration, for which the set of states where one population (or more) is extinct is not absorbing. They could express the eigenvectors of the corresponding infinitesimal generators in terms of Hahn polynomials. Many other Markov processes relevant in biological applications also admit explicit multivariate systems of polynomial eigenvectors (we refer for example to [17] where a wide class of such processes is described).

We focus here on neutral processes without mutation and immigration, which are singular for the approach of [15], and we apply our study to a much bigger class of neutral population processes, containing the birth and death processes described above, but also non-birth and death models.
Our main result is the characterization of all eigenvalues and right eigen-vectors of the transition matrix of neutral processes without mutation and immigration. To this aim, we first consider the (easier) continuous state space case in Section 2 to introduce some tools used in the sequel. Next, we construct a particular family of polynomials of two variables in Section 3, using linear algebra arguments. In Section 4, we prove that the eigenvectors of the transition matrix of neutral two-dimensional Markov processes can be decomposed as the product of “universal” polynomials (in the sense that they do not depend on the specific transition rates of the Markov chain) with functions depending only on the total population size. We then relate these eigenvectors with Dirichlet eigenvalue problems in subdomains of $Z_+^2$ of the form $\{(i, j) \in \mathbb{N}^2 : i + j \ge k\}$ for $k \ge 2$ (Section 5), where $\mathbb{N} = \{1, 2, \dots\}$.

The last section (Section 6) is devoted to the application of the previous results to the study of quasi-stationary distributions. A probability distribution $\nu$ on $Z_+^2 \setminus \{0\}$ is called quasi-stationary if it is invariant conditionally on the non-extinction of the whole population, i.e. if

$$ \mathbb{P}_{\nu}((X_1, Y_1) = (i, j) | Z_1 \neq 0) = \nu_{i,j}, \quad \forall (i, j) \in Z_+^2 \setminus \{0\}, $$

where $\mathbb{P}_{\nu}$ denotes the law of the process $(X, Y)$ with initial distribution $\nu$. This question is related to the notion of quasi-limiting distribution (also called “Yaglom limit”, in reference to Yaglom’s theorem on the same convergence for Galton-Watson processes), defined as

$$ \nu_{i,j} := \lim_{n \to +\infty} \mathbb{P}((X_n, Y_n) = (i, j) | Z_n \neq 0), \quad \forall (i, j) \in Z_+^2 \setminus \{0\}. $$

These notions are relevant in cases where extinction occurs almost surely in finite time, to describe the “stationary behaviour” of the process before extinction when the extinction time is large. This is typically the case in many population dynamics models, where ecological interactions in the population produce high mortality only when the population size is large (one speaks of density-dependent models, see e.g. [21], or [6] for discrete stochastic models).

These questions have been extensively studied in the case where the transition matrix restricted to the non-extinct states is irreducible (which is not true in our two-dimensional case). The first paper of Darroch and Seneta [7] studies the discrete-time, finite case. Several extensions of these results to continuous-time and/or infinite denumerable state spaces have then been considered in [24, 8, 11]. The case of population dynamics in dimension 1 have been studied by many authors(e.g. [5, 25, 18, 10, 16, 14, 22, 13]). More
recently, the quasi-stationary behaviour of one-dimensional diffusion models has been studied in [3]. As far as we know, the two-dimensional case has only been studied in the continuous state space (diffusion) case [4]. An extensive bibliography on quasi-stationary distributions can be found in [23]

In Subsection 6.1, we first give the quasi-limiting distribution for general finite two-dimensional Markov chains in terms of the maximal Dirichlet eigenvalues of the transition matrix in several subdomains. Finally, in Subsection 6.2, we apply our previous results to prove that coexistence in the quasi-limiting distribution is impossible for two-dimensional finite Markov chains which are close to neutrality.

The paper ends with a glossary of all the notation used in Sections 4 to 6, which may appear at different places in the paper.

## 2 Preliminary: continuous case

In this section, we consider the continuous state space case, where computations are easier, in order to introduce some of the tools needed in the discrete case.

Fix $p$ and $q$ two measurable functions from $\mathbb{R}_+$ to $\mathbb{R}_+$, and consider the system of stochastic differential equations

$$dX_t = \sqrt{2X_t p(X_t + Y_t)} dB_t^1 + X_t q(X_t + Y_t) dt \quad (1a)$$

$$dY_t = \sqrt{2Y_t p(X_t + Y_t)} dB_t^2 + Y_t q(X_t + Y_t) dt, \quad (1b)$$

where $(B^1, B^2)$ is a standard two-dimensional Brownian motion. Such SDEs are sometimes called branching diffusions, and are biologically relevant extensions of the classical Feller diffusion [3]. If $p$ and $q$ satisfy appropriate growth and regularity assumptions, the solution to this system of SDEs is defined for all positive times and can be obtained as scaling limits of two-dimensional birth and death processes (we refer to [3] for the one-dimensional case; the extension to higher dimensions is easy).

This process is neutral in the sense defined in the introduction since $Z_t = X_t + Y_t$ solves the SDE

$$dZ_t = Z_t q(Z_t) dt + \sqrt{2Z_t p(Z_t)} dB_t, \quad (2)$$

where

$$B_t = \int_0^t \sqrt{\frac{X_s}{Z_s}} dB_s^1 + \int_0^t \sqrt{\frac{Y_s}{Z_s}} dB_s^2$$
is a standard Brownian motion. Note also that $\mathbb{R}_+ \times \{0\}$, $\{0\} \times \mathbb{R}_+$ and $\{(0,0)\}$ are absorbing states as soon as there is uniqueness in law for the system (1).

For any $\varphi \in C^2(\mathbb{R}_+^2)$, the infinitesimal generator $A$ of the process $(X_t, Y_t)_{t \ge 0}$ is given by

$$ A\varphi(x, y) = \left( x \frac{\partial^2 \varphi}{\partial x^2}(x, y) + y \frac{\partial^2 \varphi}{\partial y^2}(x, y) \right) p(x + y) + \left( x \frac{\partial \varphi}{\partial x}(x, y) + y \frac{\partial \varphi}{\partial y}(x, y) \right) q(x + y). \quad (3) $$

We first observe in the following proposition that $A$ admits a symmetric measure, but only on a subset of $C^2(\mathbb{R}_+^2)$. We will use the notation $C_c^2((0, +\infty)^2)$ for the set of $C^2(\mathbb{R}_+^2)$ functions whose support is a compact subset of $(0, +\infty)^2$.

**Proposition 2.1** Assume that $1/p$ and $q/p$ belong to $L_{\text{loc}}^1((0, +\infty))$. Let us define the measure $\mu$ on $\mathbb{R}_+^2$ as

$$ \mu(dx, dy) = \frac{\exp\left(\int_1^{x+y} \frac{q(s)}{p(s)} ds\right)}{x y p(x+y)} dx dy. \quad (4) $$

Then, the restriction $\tilde{A}$ of the operator $A$ to $C_c^2((0, +\infty)^2)$ is symmetric for the canonical inner product $\langle \cdot, \cdot \rangle_\mu$ in $L^2(\mathbb{R}_+^2, \mu)$, and, hence, so is its closure in $L^2(\mathbb{R}_+^2, \mu)$.

Note that, because of the singularities in $\mu$ when $x$ or $y$ vanish, if $p \ge c > 0$ in the neighborhood of 0, any continuous function in $L^2(\mathbb{R}_+^2, \mu)$ must vanish at the boundary of $\mathbb{R}_+^2$. Therefore, $L^2(\mathbb{R}_+^2, \mu) \subset L_{\text{loc}}^{2,0}(\mathbb{R}_+^2)$, where $L_{\text{loc}}^{2,0}(\mathbb{R}_+^2)$ is defined as the closure of $C_c((0, +\infty)^2)$ in $L_{\text{loc}}^2(\mathbb{R}_+^2)$.

**Proof** For all $f, g \in C_c^2((0, +\infty)^2)$, we have (formally)

$$ \begin{aligned} \langle f, Ag \rangle_\mu = & - \int_{\mathbb{R}_+^2} p(x+y) \left( x \frac{\partial f}{\partial x}(x,y) \frac{\partial g}{\partial y}(x,y) + y \frac{\partial f}{\partial y}(x,y) \frac{\partial g}{\partial y}(x,y) \right) \mu(x,y) \, dx \, dy \\ & - \int_{\mathbb{R}_+^2} f(x,y) \left( \frac{\partial(x\mu p)}{\partial x}(x,y) \frac{\partial g}{\partial x}(x,y) + \frac{\partial(y\mu p)}{\partial y}(x,y) \frac{\partial g}{\partial y}(x,y) \right) \, dx \, dy \\ & + \int_{\mathbb{R}_+^2} f(x,y) q(x+y) \left( x \frac{\partial g}{\partial x}(x,y) + y \frac{\partial g}{\partial y}(x,y) \right) \mu(x,y) \, dx \, dy. \end{aligned} $$
Therefore, $\langle f, Ag \rangle_{\mu} = \langle Af, g \rangle_{\mu}$ if

$$xq(x + y)\mu(x, y) = \frac{\partial(x\mu p)}{\partial x}(x, y), \quad \forall x, y > 0$$

and

$$yq(x + y)\mu(x, y) = \frac{\partial(y\mu p)}{\partial y}(x, y), \quad \forall x, y > 0.$$

Conversely, these equalities can be directly checked from the formula (4), which implies that $\partial(x\mu p)/\partial x$ and $\partial(y\mu p)/\partial y$ exist in a weak sense. $\square$

Before studying the eigenvectors of A, we need the following result.

**Proposition 2.2** For all $\lambda \in \mathbb{R}$, the problem

$$(1 - x^2)h''(x) = -\lambda h(x) \qquad (5)$$

has no (weak) non-zero solution $h \in C^1([-1, 1])$ except when $\lambda = d(d-1)$ for some $d \in \mathbb{N}$. For $d = 1$, the vector space of solutions has dimension 2 and is spanned by the two polynomials $h(x) = 1$ and $h(x) = x$. For all $d \ge 2$,

$$(1 - x^2)h''(x) + d(d - 1)h(x) = 0 \qquad (6)$$

has a one-dimensional vector space of solutions in $C^1([-1, 1])$, spanned by a polynomial $H_d$ of degree $d$, which can be chosen such that the family $(H_d)_{d \ge 2}$ is an orthonormal basis of $L^2([-1, 1], \frac{dx}{1-x^2})$. In addition, $H_d$ has the same parity as $d$ (all powers of $x$ that appear are even if $d$ is even and odd if $d$ is odd). Finally, for all $d \ge 2$, one may take

$$H_d(x) = \sqrt{\frac{2d-1}{d(d-1)}} \sum_{k=1}^{d} \frac{(-d)_k (d-1)_k}{k!(k-1)!} \frac{(1+x)^k}{2^{k-1}}, \qquad (7)$$

where $(x)_n := x(x+1)\dots(x+n-1)$ is the shifted factorial.

This result is quite classical and goes back at least to Kimura [20] (up to a linear transformation), so we will omit its proof. The explicit formula for $H_d(x)$ can be obtained as a limit case of classical Gegenbauer polynomials

$$H_d^{\lambda}(x) = \sqrt{B\left(\frac{1}{2}, \lambda + \frac{1}{2}\right) \frac{2(2\lambda + 1)_d(d + \lambda)}{d!(d + 2\lambda)}} {}_2F_1\left(-d, \lambda + \frac{1}{2}; \frac{1+x}{2}\right),$$
when $\lambda \to -1/2$, where $d \ge 2$, $B$ is the Beta function and ${}_aF_b$ are the hypergeometric series. The normalizing constant in $H_d^\lambda(x)$ is not the usual one, but has been chosen such that these polynomials form an orthonormal family in $\mathbb{L}^2((-1, 1), (1-x^2)^{\lambda-1/2}dx)$.

We now introduce the change of variables

$$z = x + y \in \mathbb{R}_+ \quad \text{and} \quad w = \frac{x - y}{x + y} \in [-1, 1], \qquad (8)$$

which defines a $C^\infty$-diffeomorphism from $\mathbb{R}_+^2 \setminus \{0\}$ to $(0, +\infty) \times [-1, 1]$. Then $\mu(dx, dy) = \nu(dz, dw)$, where

$$\nu(dz, dw) := \frac{2 \exp \left( \int_1^z \frac{q(s)}{p(s)} ds \right)}{z(1 - w^2)p(z)} dz dw,$$

and

$$A = \frac{1 - w^2}{z} p(z) \frac{\partial^2}{\partial w^2} + z p(z) \frac{\partial^2}{\partial z^2} + z q(z) \frac{\partial}{\partial z}. \qquad (9)$$

Now, assume that $A$ has an eigenvector of the form $\varphi(w)\psi(z)$. If $\lambda$ is the corresponding eigenvalue, (9) yields

$$\varphi(w) (zp(z)\psi''(z) + zq(z)\psi'(z) - \lambda\psi(z)) + (1 - w^2)\varphi''(w) \frac{p(z)}{z} \psi(z) = 0.$$

Hence, we have the following result.

**Proposition 2.3** All functions on $\mathbb{R}_+^2$ of the form

$$H_d\left(\frac{x-y}{x+y}\right)\psi(x+y), \qquad (10)$$

where $d \ge 0$, $H_1(x) = x$, $H_0(x) = 1$, $H_k$, $k \ge 2$ are as in Proposition 2.2, where $\psi$ satisfies

$$zp(z)\psi''(z) + zq(z)\psi'(z) - d(d-1)\frac{p(z)}{z}\psi(z) = \lambda\psi(z), \quad \forall z \ge 0, \qquad (11)$$

for some $\lambda \in \mathbb{R}$, are eigenfunctions of $A$ for the eigenvalue $\lambda$.
Now, we proceed to prove that there are no other eigenvectors of A (in the sense of Theorem 2.4 below). Such a result seems natural as (11) can be written as the Sturm-Liouville problem

$$ \frac{d}{dz} \left( \psi'(z) \exp \int_{1}^{z} \frac{q(s)}{p(s)} ds \right) - d(d-1) \frac{\exp \int_{1}^{z} \frac{q(s)}{p(s)} ds}{z^2} \psi(z) = -\lambda \frac{\exp \int_{1}^{z} \frac{q(s)}{p(s)} ds}{z p(z)} \psi(z). $$

Here again, the usual integrability conditions are not satisfied. More precisely, if $p(z) \ge c > 0$ for $z$ in the neighborhood of 0 and $q$ is bounded, using the terminology of [26], for all $d \ge 0$, the problem (10) is a singular self-adjoint boundary value problem on $(0, \infty)$, where the endpoint 0 is LP (limit point) singular (see [26, Thm. 7.4.1]). In this case very little is known about the existence of an orthonormal basis of eigenvectors in $L^2((0, \infty), \tilde{\nu})$, where

$$ \tilde{\nu}(z) = \frac{2 \exp \left( \int_{1}^{z} \frac{q(s)}{p(s)} ds \right)}{z p(z)} \quad (12) $$

(the spectrum might even be continuous, see [26, Thm. 10.12.1(8)]).

For this reason, we state our next result on the operator A with a restricted domain corresponding to the case where the diffusion is reflected in the set

$$ D := \{(x, y) \in \mathbb{R}_+^2 : a \le x + y \le b\}, $$

where $0 < a < b < \infty$. For all initial condition $(X_0, Y_0) = (x_0, y_0) \in D$, we consider the process $(X_t, Y_t, k_t^{(a)}, k_t^{(b)})_{t \ge 0}$ such that $k_0^{(a)} = k_0^{(b)} = 0$, $k^{(a)}$ and $k^{(b)}$ are nondecreasing processes, $(X_t, Y_t) \in D$ for all $t \ge 0$,

$$ k_t^{(\alpha)} = \int_0^t 1_{\{X_s+Y_s=a\}} dk_s^{(\alpha)}, \quad \forall t \ge 0, \quad \alpha = a, b $$

$$ dX_t = \sqrt{2X_t p(X_t + Y_t)} dB_t^1 + X_t q(X_t + Y_t) dt - \sqrt{2}^{-1} dk_t^{(b)} + \sqrt{2}^{-1} dk_t^{(a)} $$

$$ dY_t = \sqrt{2Y_t p(X_t + Y_t)} dB_t^2 + Y_t q(X_t + Y_t) dt - \sqrt{2}^{-1} dk_t^{(b)} + \sqrt{2}^{-1} dk_t^{(a)}. $$

Then $Z_t = X_t + Y_t$ is the solution of (2) reflected at $a$ and $b$ with local time $k_t^{(a)}$ at $a$ (resp. $k_t^{(b)}$ at $b$). Therefore, $(X_t, Y_t)$ is also neutral in the sense of the introduction.

The corresponding infinitesimal generator is defined by (3) with domain the set of $\varphi(x, y) \in C^1(D) \cap C^2(\text{int}(D))$, where $\text{int}(D)$ denotes the interior of $D$, such that

$$ \frac{\partial \varphi}{\partial x}(x, y) + \frac{\partial \varphi}{\partial y}(x, y) = 0, \quad \forall (x, y) \in D \text{ s.t. } x + y = a \text{ or } b. $$
**Theorem 2.4** For $0 < a < b < \infty$, assume that $p \ge c > 0$ on $[a, b]$ and $q/p$ belong to $L^1([a, b])$.

(a) There exists a denumerable orthonormal basis of $L^2(D, \mu)$ of eigenvectors of $A$ of the form (10), where $d \ge 2$ and $\psi$ solves (11) on $(a, b)$ and satisfies $\psi'(a) = \psi'(b) = 0$. Moreover, any eigenvector of $A$ in $L^2(D, \mu)$ is a linear combination of eigenvectors of the form (10), all corresponding to the same eigenvalue.

(b) There exists a family of right eigenvectors of $A$ of the form

$$ \{1\} \cup \left\{ \frac{x}{x+y} \psi_k(x+y), \frac{y}{x+y} \psi_k(x+y) \right\}_{k \ge 1} \cup \bigcup_{d \ge 2} \left\{ H_d \left( \frac{x-y}{x+y} \right) \psi_k^{(d)}(x+y) \right\}_{k \ge 1}, \quad (13) $$

which is a basis of the vector space

$$ V := \left\{ f \in L^2(D, \text{Leb}) : \exists f_1, f_2 \in L^2([a, b], \text{Leb}) \text{ and } f_3 \in L^2(D, \mu) \text{ s.t.} \\ f(x, y) = \frac{x}{x+y} f_1(x+y) + \frac{y}{x+y} f_2(x+y) + f_3(x, y) \right\}, \quad (14) $$

where Leb denotes Lebesgue's measure. More precisely, for all $f \in V$, the functions $f_1, f_2, f_3$ in (14) are unique and there exists unique sequences $\{\alpha_k\}_{k \ge 1}$, $\{\beta_k\}_{k \ge 1}$, and $\{\gamma_{dk}\}_{d \ge 2, k \ge 1}$ such that

$$ f(x, y) = \sum_{k \ge 1} \alpha_k \frac{x}{x+y} \psi_k(x+y) + \sum_{k \ge 1} \beta_k \frac{y}{x+y} \psi_k(x+y) \\ + \sum_{d \ge 2} \sum_{k \ge 1} \gamma_{dk} H_d \left( \frac{x-y}{x+y} \right) \psi_k^{(d)}(x+y), \quad (15) $$

where the series $\sum_k \alpha_k \psi_k$ and $\sum_k \beta_k \psi_k$ both converge for $\| \cdot \|_{\tilde{\nu}}$ and $\sum_{d,k} \gamma_{dk} H_d(\frac{x-y}{x+y}) \psi_k^{(d)}(x+y)$ converges for $\| \cdot \|_{\mu}$.

Point (b) says that the eigenvectors of the form (10), although not orthogonal in some Hilbert space, allow one to recover a bigger class of functions than in Point (a). The vector space $V$ is not equal to $L^2(D, \text{Leb})$, but the following result shows that it is much bigger than $L^2(D, \mu)$.

**Proposition 2.5** The vector space $V$ of Theorem 2.4 contains $H^1(D, \text{Leb})$.
**Proof** To prove this result, it is more convenient to consider the variables $(z, w)$ as in (8) instead of $(x, y)$. The vector space $V$ then becomes the set of $g \in \mathbb{L}^2([a, b] \times [-1, 1], \text{Leb})$ such that

$$g(z, w) = \frac{1+w}{2}g_1(z) + \frac{1-w}{2}g_2(z) + g_3(z, w) \quad (16)$$

for some $g_1, g_2 \in \mathbb{L}^2([a, b], \text{Leb})$ and $g_3 \in \mathbb{L}^2([a, b] \times [-1, 1], \nu) = \mathbb{L}^2([a, b] \times [-1, 1], (1-w^2)^{-1}dz dw)$.

We first introduce the following notion of trace: we say that a function $g \in \mathbb{L}^2([a, b] \times [-1, 1], \text{Leb})$ admits the function $\bar{g} \in \mathbb{L}^2([a, b], \text{Leb})$ as a trace at $w = 1$, or $w = -1$ respectively, if

$$g(z, w) - \bar{g}(z) \in \mathbb{L}^2([a, b] \times [0, 1], (1-w^2)^{-1}dz dw),$$

or

$$g(z, w) - \bar{g}(z) \in \mathbb{L}^2([a, b] \times [-1, 0], (1-w^2)^{-1}dz dw)$$

respectively.

Our first claim is that any $g \in \mathbb{L}^2([a, b] \times [-1, 1], \text{Leb})$ which admits traces $g_1$ and $g_2$ at $w = 1$ and $w = -1$ respectively, belongs to $V$, and these traces are exactly the functions $g_1$ and $g_2$ in (16). To see that, we only have to check that $g_3 \in \mathbb{L}^2([a, b] \times [0, 1], (1-w^2)^{-1}dz dw)$, and the same result on $[a, b] \times [-1, 0]$ will follow by symmetry:

$$\int_a^b dz \int_0^1 dw \frac{\left(g - \frac{1+w}{2}g_1 - \frac{1-w}{2}g_2\right)^2}{1-w^2} \\ \le 2 \int_a^b dz \int_0^1 dw \frac{(g-g_1)^2}{1-w} + 2 \int_a^b dz \int_0^1 dw (g_1+g_2)^2 \frac{1-w}{4(1+w)} < +\infty.$$

Second, we claim that any $g \in H^1([a, b] \times [-1, 1], \text{Leb})$ admits traces at $w = 1$ and $w = -1$ as defined above. Assume first that $g \in C^1([a, b] \times [0, 1])$. Then, using the Cauchy-Schwartz inequality,

$$\int_0^1 dw \int_a^b dz \frac{(g(z, w) - g(z, 1))^2}{1 - w^2} = \int_0^1 dw \int_a^b dz \frac{\left(\int_w^1 \nabla_w g(z, x) dx\right)^2}{1 - w^2} \\ \le \int_0^1 dw \int_a^b dz \int_w^1 dx |\nabla_w g(z, x)|^2 \le \|\nabla_w g\|_{\text{Leb}}^2.$$
Since in addition $\|g(\cdot, 1)\|_{L^2(\text{Leb})} \le 4\|g\|_{H^1(\text{Leb})}$ by classical trace results (cf. e.g. [2, p. 196]), the function $g \mapsto (g(\cdot, 1), g - g(\cdot, 1))$ extends by density to a linear operator $\psi$ continuous from $H^1([a, b] \times [0, 1], \text{Leb})$ to $L^2([a, b], \text{Leb}) \times L^2([a, b] \times [0, 1], (1 - w^2)^{-1} dz dw)$. Since obviously $\psi_1(g) + \psi_2(g) = g$, the claim is proved and the proof of Proposition 2.5 is completed. $\square$

**Proof of Theorem 2.4** An eigenvector of A of the form $H_d(w)\psi(z)$ satisfies the Neumann boundary condition in D iff $\psi'(a) = \psi'(b) = 0$. The problem (11) with this boundary condition is a regular Sturm-Liouville problem with the weight $\tilde{\nu}$ defined in (12). Therefore (cf. e.g. [26, Thm. 4.6.2]), for all $d \ge 0$, there exists an orthonormal basis $\{\psi_k^{(d)}\}_{k \ge 1}$ of $L^2([a, b], \tilde{\nu})$ composed of solutions to (11) on $(a, b)$ with Neumann boundary conditions. All the corresponding eigenvalues are real, simple and the corresponding spectrum has no accumulation point.

Now, we claim that

$$ \mathcal{F} := \bigcup_{d \ge 2} \left\{ H_d \left( \frac{x-y}{x+y} \right) \psi_k^{(d)} (x+y) \right\}_{k \ge 1} $$

forms an orthonormal basis of $L^2(D, \mu)$. The orthonormal property follows from the fact that, if $\varphi(x, y) = H_d(\frac{x-y}{x+y})\psi_k^{(d)}(x+y)$ and $\varphi'(x, y) = H_{d'}(\frac{x-y}{x+y})\psi_{k'}^{(d')}(x+y)$ for $d, d' \ge 2$ and $k, k' \ge 1$,

$$ \begin{aligned} \langle \varphi, \varphi' \rangle_{\mu} &= \int_{-1}^{1} \int_{a}^{b} H_{d}(w) H_{d'}(w) \psi_{k}^{(d)}(z) \psi_{k'}^{(d')}(z) d\nu(z, w) \\ &= \langle \psi_{k}^{(d)}, \psi_{k'}^{(d')} \rangle_{\tilde{\nu}} \langle H_{d}, H_{d'} \rangle_{(1-w^2)^{-1} dw}. \end{aligned} $$

To prove that $\mathcal{F}$ is a basis of $L^2(D, \mu)$, assume that $f \in L^2(D, \mu)$ satisfies

$$ \iint_D f(x, y) H_d \left( \frac{x-y}{x+y} \right) \psi_k^{(d)} (x+y) d\mu(x, y) = 0, \quad \forall d \ge 2, k \ge 1. $$

By the Cauchy-Schwartz inequality, for all $d \ge 2$, the function

$$ \tilde{f}_d(z) := \int_{-1}^{1} f \left( \frac{z(1+w)}{2}, \frac{z(1-w)}{2} \right) \frac{H_d(w)}{1-w^2} dw $$
belongs to $\mathbb{L}^2([a, b], \tilde{\nu})$. In addition,

$$ \langle \tilde{f}_d, \psi_k^{(d)} \rangle_{\tilde{\nu}} = 0, \quad \forall k \geq 1. $$

Therefore, $\tilde{f}_d(z) = 0$ for all $d \geq 2$, for Lebesgue-almost every $z \geq 0$. By Fubini's theorem, $w \mapsto f(\frac{z(1+w)}{2}, \frac{z(1-w)}{2})$ belongs to $\mathbb{L}^2([-1, 1], (1-x^2)^{-1}dx)$ for almost every $z \geq 0$. Hence we deduce from Proposition 2.2 that this function is 0 for almost every $z \geq 0$. Hence $f = 0$.

Thus $\mathcal{F}$ is an orthonormal basis of $\mathbb{L}^2(D, \mu)$ composed of eigenvectors of $A$. It is then classical to deduce that $A$ admits no other eigenvector in this space, in the sense of point (a).

For point (b), let us first prove that the decomposition

$$ f = \frac{x}{x+y}f_1(x+y) + \frac{y}{x+y}f_2(x+y) + f_3(x,y) $$

is unique for $f \in V$, with $f_1, f_2 \in \mathbb{L}^2([a, b], \text{Leb})$ and $f_3 \in \mathbb{L}^2(D, \mu)$. We only need to prove that this equality for $f = 0$ implies $f_1 = f_2 = f_3 = 0$.

Since $f_3 \in \mathbb{L}^2(D, \mu)$, we have

$$ \int_{0}^{\varepsilon} dx \int_{a-x}^{b-x} dy f_3^2(x, y) \leq \varepsilon \int_{0}^{\varepsilon} dx \int_{a-x}^{b-x} dy \frac{f_3^2(x, y)}{x} = o(\varepsilon) $$

as $\varepsilon \to 0$. Therefore,

$$ \begin{aligned} \int_{0}^{\varepsilon} dx \int_{a-x}^{b-x} dy \frac{y^2}{(x+y)^2} f_2^2 &= \int_{0}^{\varepsilon} dx \int_{a-x}^{b-x} dy \left( \frac{x}{x+y} f_1^2 + f_3 \right)^2 \\ &\leq 2 \int_{0}^{\varepsilon} dx \int_{a-x}^{b-x} dy \frac{x^2}{(x+y)^2} f_1^2 + o(\varepsilon) = o(\varepsilon). \end{aligned} $$

This implies that $\int_a^b f_2^2(z)dz = 0$, i.e. $f_2 = 0$. Similarly, $f_1 = 0$ and thus $f_3 = 0$.

Since $\mathbb{L}^2([a, b], \text{Leb}) = \mathbb{L}^2([a, b], \tilde{\nu})$, the result then follows from the decomposition of $f_1$ and $f_2$ (resp. $f_3$) in the orthonormal basis $\{\psi_k^{(1)}\}_{k \geq 1}$ of $\mathbb{L}^2([a, b], \tilde{\nu})$ (resp. $\{H_d(\frac{x-y}{x+y})\psi_k^{(d)}(x+y)\}_{d \geq 2, k \geq 1}$ of $\mathbb{L}^2(D, \mu)$). $\square$

To motivate the calculations of the next section, let us finally observe that, for all $\varphi \in C^2(\mathbb{R}_+^2)$,

$$ A\varphi(x, y) = (\tilde{T}\varphi)(x, y)p(x + y) + (\tilde{L}\varphi)(x, y)q(x + y), $$
where

$$ \tilde{L} = x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} $$

and

$$ \tilde{T} = x \frac{\partial^2}{\partial x^2} + y \frac{\partial^2}{\partial y^2}, \qquad (17) $$

and that $\tilde{T}\tilde{L} = \tilde{L}\tilde{T}$.

### 3 On a family of bivariate polynomials

The goal of this section is to prove the existence of a family of polynomials in $\mathbb{R}$ of two variables $X$ and $Y$, satisfying the family of relations

$$ XP(X+1, Y) + YP(X, Y+1) = (X + Y + d)P(X, Y) \qquad (18a) $$

$$ XP(X-1, Y) + YP(X, Y-1) = (X + Y - d)P(X, Y) \qquad (18b) $$

for an integer $d \ge 0$.

Before stating the main result of the section, let us recall some notation. $\mathbb{R}[X]$ is the set of polynomials on $\mathbb{R}$ with a single variable $X$ and $\mathbb{R}[X, Y]$ the set of real polynomials with two variables $X$ and $Y$. The degree $\deg(P)$ of a polynomial $P \in \mathbb{R}[X, Y]$ is defined as the maximal total degree of each monomial of $P$. We define

$$ \mathcal{P}_d = \{P \in \mathbb{R}[X, Y] : \deg(P) \le d\}. $$

For all $P \in \mathbb{R}[X, Y]$, we may write

$$ P(X, Y) = \sum_{i,j \ge 0} a_{i,j} X^i Y^j, $$

where only finitely many of the $a_{i,j}$ are nonzero. The real number $a_{i,j}$ will be called the $(i, j)$-coefficient of $P$.

For any $P \in \mathbb{R}[X, Y]$ and for any $d \ge 0$, we denote by $[P]_d$ the sum of all monomials of $P$ of degree $d$:

$$ [P]_d(X, Y) = \sum_{i=0}^{d} a_{i,d-i} X^i Y^{d-i}. $$
In particular, $[P]_d$ is homogeneous of degree $d$ and $P = \sum_{i=1}^{\infty}[P]_i$.

We denote by $\Delta_i$ the first-order symmetrized discrete derivative with respect to the $i$-th variable:

$$ \forall P \in \mathbb{R}[X, Y], \quad \Delta_1 P(X, Y) = \frac{P(X+1, Y) - P(X-1, Y)}{2} \\ \text{and} \quad \Delta_2 P(X, Y) = \frac{P(X, Y+1) - P(X, Y-1)}{2}, $$

and by $\Delta_i^2$ the symmetrized second-order discrete derivative with respect to the $i$-th variable:

$$ \forall P \in \mathbb{R}[X, Y], \quad \Delta_1^2 P(X, Y) = \frac{P(X+1, Y) + P(X-1, Y) - 2P(X, Y)}{2} \\ \text{and} \quad \Delta_2^2 P(X, Y) = \frac{P(X, Y+1) + P(X, Y-1) - 2P(X, Y)}{2}. $$

Note that the superscript in the notation $\Delta_i^2$ does not correspond to the composition of the operator $\Delta_i$ with itself.

Finally, we define the linear operators on $\mathbb{R}[X, Y]$

$$ L = X\Delta_1 + Y\Delta_2 \quad \text{and} \quad T = X\Delta_1^2 + Y\Delta_2^2. $$

Then, adding and subtracting the equations (18a) and (18b), the system (18) is equivalent to

$$ LP = dP \qquad (19) $$

$$ TP = 0. \qquad (20) $$

We are going to prove the following result

**Theorem 3.1** For $d = 1$, the system (18) has a two-dimensional vector space of solutions in $\mathbb{R}[X, Y]$, spanned by the two polynomials $P_1^{(1)} := X$ and $P_1^{(2)} := Y$.

For any $d \in \{0, 2, 3, \dots\}$, the system (18) has a one-dimensional vector space of solutions. All nonzero solutions are of degree $d$. For $d = 0$, this is the vector space of constants, spanned by $P_0 := 1$. When $d \ge 2$, we denote by $P_d$ the unique solution to (18) with $(d-1, 1)$-coefficient equal to $-2H_d'(1) = (-1)^d 2\sqrt{d(d-1)(2d-1)}$, where $H_d$ is defined in Proposition 2.2.
It can be checked that the first polynomials are

$$
\begin{align*}
P_0 &= 1, \quad P_1^{(1)} = X, \quad P_1^{(2)} = Y \\
P_2 &= 2\sqrt{6} XY, \\
P_3 &= -2\sqrt{30} XY(X-Y), \\
P_4 &= 4\sqrt{21} XY(X^2 - 3XY + Y^2 + 1), \\
P_5 &= -6\sqrt{20} XY(X-Y)(X^2 - 5XY + Y^2 + 5).
\end{align*}
$$

Before proving this result, let us give some properties of the polynomials $P_d$, proved after the proof of Theorem 3.1.

**Proposition 3.2** The polynomials $P_d$, $d \ge 2$, defined in Theorem 3.1 satisfy
the following properties:

(a) For all $d \ge 2$, $[P_d]_d(X, Y) = (X+Y)^d H_d(\frac{X-Y}{X+Y})$, where $H_d$ is defined in Proposition 2.2.

(b) $[P]_{d-2k-1} = 0$ for all $0 \le k < d/2$.

(c) For all $d \ge 2$, $P_d$ is divisible by $XY$. For $d$ odd, $P_d(X, Y)$ is divisible by $XY(X - Y)$.

(d) for all $d \ge 2$, $P_d(Y, X) = P_d(-X, -Y) = (-1)^d P_d(X, Y)$.

(e) $P_d(i, j) = 0$ if $i, j \in \mathbb{Z}$, $ij \ge 0$ and $0 \le |i| + |j| \le d - 1$.

(f) For all $d \ge 0$, the matrix $(P_i(j, d-j))_{0 \le i, j \le d}$ is invertible, where $P_1 = P_1^{(1)}$.
In particular, $(P_d(j, d-j))_{0 \le j \le d} \ne 0$.

(g) For all $d \ge 3$, $P_d(j, d-j)P_d(j+1, d-j-1) < 0$ if $1 \le j \le d-2$.

(h) For all $d, d', k \ge 2$,

$$
\sum_{i=1}^{k-1} \frac{P_d(i, k-i) P_{d'}(i, k-i)}{i(k-i)} = 2 \binom{k+d-1}{2d-1} \delta_{dd'}, \quad (21)
$$

where $\delta_{ij}$ is the Kronecker symbol and by convention $\binom{i}{j} = 0$ if $j < 0$
or $j > i$.
**Proposition 3.3** For all $d \ge 2$, the polynomial $P_d$ is given by the following formula:

$$P_d(X, Y) = C_d(-X - Y)_d \sum_{k=1}^{d} \frac{(-d)_k (d-1)k}{(k-1)!k!} \frac{(-X)_k}{(-X-Y)_k}, \quad (22)$$

where

$$C_d = (-1)^{d+1} 2 \sqrt{\frac{2d-1}{d(d-1)}}.$$

In order to prove Theorem 3.1, we need the following lemma.

**Lemma 3.4**

(a) We have

$$TL = LT + T, \quad (23)$$

(b) Define for all $d \ge 0$

$$\mathcal{D}_d = \{P \in \mathcal{P}_d : T(P) = 0\}.$$

Then $\dim(\mathcal{D}_0) = 1$ and $\dim(\mathcal{D}_d) = d + 2$ for all $d \ge 1$.

**Proof of Lemma 3.4** To prove (23), it would be enough to expand both sides of the equation. We prefer to give a proof based on differential calculations, because it is related to the method used in the rest of the proof. First, let $I = \{1, 3, 5, \dots\}$ and $J = \{2, 4, 6, \dots\}$ be the sets of odd and even integers, respectively. Using the fact that, for $Q \in \mathbb{R}[X]$.

$$Q(X + 1) = \sum_{n \ge 0} \frac{Q^{(n)}(X)}{n!} \quad \text{and} \quad Q(X - 1) = \sum_{n \ge 0} \frac{(-1)^n Q^{(n)}(X)}{n!},$$

where $Q^{(n)}$ denotes the $n$-th derivative of $Q$, one has

$$T = \sum_{p \in J} \frac{1}{p!} \left( X \frac{\partial^p}{\partial X^p} + Y \frac{\partial^p}{\partial Y^p} \right), \quad (24)$$

$$L = \sum_{q \in I} \frac{1}{q!} \left( X \frac{\partial^q}{\partial X^q} + Y \frac{\partial^q}{\partial Y^q} \right). \quad (25)$$
Since for all $p, q \in \mathbb{N}$

$$X \frac{\partial^p}{\partial X^p} \left( X \frac{\partial^q}{\partial X^q} \right) = X^2 \frac{\partial^{p+q}}{\partial X^{p+q}} + pX \frac{\partial^{p+q-1}}{\partial X^{p+q-1}},$$

one easily checks that

$$TL - LT = \sum_{p \in J, q \in I} \frac{1}{p!} \frac{1}{q!} (p-q) \left( X \frac{\partial^{p+q-1}}{\partial X^{p+q-1}} + Y \frac{\partial^{p+q-1}}{\partial Y^{p+q-1}} \right).$$

Now, for all $n \in \mathbb{N}$,

$$\begin{align*}
\sum_{p \in J, q \in I, p+q=2n+1} \frac{1}{p!} \frac{1}{q!} (p-q) &= \sum_{p \in I, q \in I, p+q=2n} \frac{1}{p!} \frac{1}{q!} - \sum_{p \in J, q \in J \cup \{0\}, p+q=2n} \frac{1}{p!} \frac{1}{q!} \\
&= - \sum_{0 \le p \le 2n} \frac{1}{p!} \frac{1}{(2n-p)!} (-1)^p + \frac{1}{(2n)!} \\
&= ((1-1)^{2n} + 1) \frac{1}{(2n)!} = \frac{1}{(2n)!}.
\end{align*}$$

This completes the proof of Lemma 3.4 (a).

To prove (b), let $\mathcal{H}_p$ be the subspace of $\mathcal{P}_d$ composed of all homogeneous polynomials of degree $d$. Recall the definition (17) of the operator $\tilde{T}$ on $\mathbb{R}[X, Y]$, and observe that, for all $d \ge 1$, $\tilde{T}$ is a linear map from $\mathcal{H}_d$ to $\mathcal{H}_{d-1}$. Now, the family $\{(X-Y)^k(X+Y)^{d-k}, 0 \le k \le d\}$ forms a basis of $\mathcal{H}_d$. Hence, any $P \in \mathcal{H}_d$ can be written in the form $P(X, Y) = (X+Y)^d h(\frac{X-Y}{X+Y})$, where $h \in \mathbb{R}[X]$ has degree $d$. With this notation, it can be checked that

$$\begin{align*}
\frac{\partial^2 P}{\partial X^2}(X, Y) &= d(d-1)(X+Y)^{d-2}h(W) \\
&\quad + 4(d-1)Y(X+Y)^{d-3}h'(W) + 4Y^2(X+Y)^{d-4}h''(W),
\end{align*}$$

where $W = (X - Y)/(X + Y)$, and similarly for the second variable. This yields

$$\tilde{T}P = (X + Y)^{d-1} \left( d(d-1)h(W) + 4 \frac{XY}{(X+Y)^2} h''(W) \right).$$

Using the relation $4XY/(X+Y)^2 = 1 - W^2$, we finally obtain that $P \in \text{Ker}(\tilde{T}) \cap \mathcal{H}_d$ if and only if $h$ solves (6). By Proposition 2.2, for all $d \ne 1$, this
equation has a unique (up to a multiplicative constant) polynomial solution, which has degree $d$. Since $\dim(\mathcal{H}_d) = d+1$, we deduce that $\tilde{T}: \mathcal{H}_d \to \mathcal{H}_{d-1}$ is surjective for all $d \ge 2$. If $d=1$, $\ker(\tilde{T}) \cap \mathcal{H}_1 = \mathcal{H}_1$ which has dimension 2.

Now, let $P = [P]_1 + \dots + [P]_d \in \mathcal{P}_d$ and observe that any $k$-th order derivative of $[P]_m$ belongs to $\mathcal{H}_{m-k}$ if $k \le m$. Therefore, by (24), the equation $TP = 0$ is equivalent to the fact that, for all $0 \le n \le d-1$,

$$[TP]_n = \sum_{p \ge 1} \frac{1}{(2p)!} \left( X \frac{\partial^{2p}}{\partial X^{2p}} + Y \frac{\partial^{2p}}{\partial Y^{2p}} \right) [P]_{n+2p-1} = 0,$$

or, equivalently,

$$\tilde{T}[P]_{n+1} = -2 \sum_{p \ge 2} \frac{1}{(2p)!} \left( X \frac{\partial^{2p}}{\partial X^{2p}} + Y \frac{\partial^{2p}}{\partial Y^{2p}} \right) [P]_{n+2p-1}. \quad (26)$$

If $n \ge 1$ and $[P]_{n+3}, [P]_{n+5}, \dots$ are given, there is a one-dimensional affine space of solution for this equation. If $n=0$, (26) is automatically satisfied, since both sides are 0. Therefore, choosing recursively $[P]_d, [P]_{d-1}, \dots, [P]_2$ and setting any value to $[P]_1$ and $[P]_0$, the result on the dimension of $\mathcal{D}_d$ easily follows. $\square$

**Proof of Theorem 3.1** Fix $d \ge 0$. We claim that, as a linear operator on $\mathcal{P}_d$, $L$ is diagonalizable and its spectrum $\text{Sp}_{\mathcal{P}_d}(L) = \{0, 1, \dots, d\}$. To see this, fix $\lambda \in \text{Sp}_{\mathcal{P}_d}(L)$ and $P$ an eigenvector for this eigenvalue, with degree $p$. Writing as in the proof of Lemma 3.4 $P = [P]_p + \dots + [P]_0$, the equation $LP = \lambda P$ is equivalent to the fact that, for $0 \le n \le p$

$$\lambda[P]_n = \sum_{q \ge 0} \frac{1}{(2q+1)!} \left( X \frac{\partial^{2q+1}}{\partial X^{2q+1}} + Y \frac{\partial^{2q+1}}{\partial Y^{2q+1}} \right) [P]_{n+2q}. \quad (27)$$

Now, for any $Q \in \mathcal{H}_k$, one has

$$\tilde{L}Q = X \frac{\partial Q}{\partial X} + Y \frac{\partial Q}{\partial Y} = kQ. \quad (28)$$

Therefore, (27) for $n=p$ imposes $\lambda=p$, and for $n=p-1$, $[P]_{p-1} = 0$. Moreover, for $0 \le n \le p-2$, (27) is equivalent to

$$(p-n)[P]_n = \sum_{q \ge 1} \frac{1}{(2q+1)!} \left( X \frac{\partial^{2q+1}}{\partial X^{2q+1}} + Y \frac{\partial^{2q+1}}{\partial Y^{2q+1}} \right) [P]_{n+2q}, \quad (29)$$
which allows one to compute recursively $[P]_{p-2}, \dots, [P]_0$ given any $[P]_p \in \mathcal{H}_p$ and $[P]_{p-1} = 0$. Since $\dim(\mathcal{H}_p) = p+1$, the eigenspace corresponding to the eigenvalue $p$ of $L$ has dimension $p+1$.

Now, it follows from (23) that $L$ is a linear operator from $\mathcal{D}_d$ to $\mathcal{D}_d$. Since $L$ is diagonalizable in $\mathcal{P}_d$, it is also diagonalizable on any stable subspace, and $\mathrm{Sp}_{\mathcal{D}_d}(L) \subset \{0, \dots, d\}$.

Let $p \in \{2, \dots, d\}$ and assume that there exists $P \in \mathcal{D}_d \setminus \{0\}$ satisfying $LP = pP$. Again, $\deg(P) = p$ necessarily. Writing $P = [P]_p + \dots + [P]_0$ again, since $T(P) = 0$, we have

$$ \tilde{T}[P]_p = X \frac{\partial^2[P]_p}{\partial X^2} + Y \frac{\partial^2[P]_p}{\partial Y^2} = 0, $$

which has a one-dimensional vector space of solutions in $\mathcal{H}_p$. Once $[P]_p$ is fixed and since we have $[P]_{p-1} = 0$, (29) can be used recursively to compute $[P]_{p-2}, \dots, [P]_0$. In conclusion, the eigenspace of $L$ in $\mathcal{D}_d$ for the eigenvalue $p$ is either of dimension 1 or 0. Now, $L$ is diagonalizable in $\mathcal{D}_d$. Since $\dim(\mathcal{D}_d) = d+2$ and $\dim(\mathcal{H}_0) + \dim(\mathcal{H}_1) = 3$, the only possibility is that $\mathrm{Sp}_{\mathcal{D}_d}(L) = \{0, 1, \dots, d\}$ and that each eigenvalue $p \neq 1$ has a one-dimensional vector space of solutions, and the eigenvalue 1 has a two-dimensional vector space of solutions.

This easily implies Theorem 3.1, except for the expression of $-2H'_d(1)$. This can be easily obtained from (7) and the relation $H'_d(1) = (-1)^{d-1}H'_d(-1)$ which follows from the parity property of $H_d$ stated in Proposition 2.2. $\square$

**Proof of Proposition 3.2** Recall from the proof of Lemma 3.4 that a nonzero solution $P$ of (18) for $d \geq 2$ satisfies

$$ [P]_d(X, Y) = \sum_{n=0}^{d} b_n (X - Y)^n (X + Y)^{d-n}, $$

where $h(x) = \sum_{n=0}^{d} b_n x^n$ is a polynomial solution of (6). Therefore, the $(d-1, 1)$-coefficient $a_{d-1, 1}$ of $P$ is given by

$$ a_{d-1, 1} = \sum_{n=0}^{d} b_n (-n + d - n) = dh(1) - 2h'(1) = -2h'(1), $$
and Point (a) then follows from Proposition 2.2 and the value of the $(d-1, 1)$ coefficient of $P_d$.

Observe that any polynomial solution of (6) with $d \ge 2$ is divisible by $(X-1)(X+1)$. As a consequence of the previous construction, any polynomial $P$ such that (18) holds satisfies that $XY$ divides $[P]_d$. Note also that $[P]_{d-1} = 0$, which implies (b) by (29). Moreover, (29) also implies by induction that $XY$ divides $[P]_{d-2k}$ for all $0 \le k \le d/2$, which yields the first part of (c).

By Proposition 2.2, (d) is true for $[P_d]_d$ and of course for $[P_d]_{d-1} = 0$. Now, assume that $P \in \mathbb{R}[X, Y]$ satisfies $P(Y, X) = P(-X, -Y) = (-1)^\alpha P(X, Y)$. Then it can be easily checked that, for all $k \ge 1$,

$$Q(X, Y) := X \frac{\partial^k P}{\partial X^k}(X, Y) + Y \frac{\partial^k P}{\partial Y^k}(X, Y)$$

satisfies

$$Q(Y, X) = (-1)^\alpha Q(X, Y) \quad \text{and} \quad Q(-X, -Y) = (-1)^{\alpha+k+1} Q(X, Y).$$

Therefore, (d) easily follows from (29) by induction.

Now, fix $d$ odd. By Proposition 2.2, $H_d$ is odd and thus $H_d(X)$ is divisible by $X$. This implies that $[P]_d$ is divisible by $X-Y$. Moreover, it follows from (d) that the polynomial $[P]_d/(X-Y)$ is symmetric. Now, let $Q(X, Y) = X^n Y^m + X^m Y^n$ for some $n, m \ge 0$, and fix $k \ge 0$. Since

$$X \frac{\partial^k Q}{\partial X^k} - Y \frac{\partial^k Q}{\partial Y^k} = n(n-1) \dots (n-k+1)(X^{n-k+1} Y^m - X^m Y^{n-k+1}) \\ + m(m-1) \dots (m-k+1)(X^{m-k+1} Y^n - X^n Y^{m-k+1}),$$

the polynomial $X-Y$ divides $X \frac{\partial^k Q}{\partial X^k} - Y \frac{\partial^k Q}{\partial Y^k}$. Since this holds for any $n, m \ge 0$, the same is true for all $Q$ such that $Q(X, Y) = Q(Y, X)$. Now, any polynomial of the form $P(X, Y) = (X-Y)Q(X, Y)$ with $Q$ symmetric satisfies

$$X \frac{\partial^k P}{\partial X^k}(X, Y) + Y \frac{\partial^k P}{\partial Y^k}(X, Y) = kX \frac{\partial^{k-1} Q}{\partial X^{k-1}}(X, Y) - kY \frac{\partial^{k-1} Q}{\partial Y^{k-1}}(X, Y) \\ + (X-Y) \left( X \frac{\partial^k Q}{\partial X^k}(X, Y) + Y \frac{\partial^k Q}{\partial Y^k}(X, Y) \right).$$

In particular, $X-Y$ divides $X \frac{\partial^k P}{\partial X^k} + Y \frac{\partial^k P}{\partial Y^k}$. Therefore, the fact that $[P]_i$ is divisible by $X-Y$ for $i < d$ follows from (29) by induction. This ends the proof of (c).
As a consequence of (c), $P_d(i, 0) = P_d(0, j) = 0$ for any $i, j \in \mathbb{Z}$ for $d \ge 2$. Applying (18b) for $(X, Y) = (d-1, 1)$ yields $P_d(d-2, 1) = 0$. By induction, applying (18b) for $(X, Y) = (d-k, 1)$ implies $P_d(d-k-1, 1) = 0$ for all $k \in \{1, \dots, d-2\}$. Similarly, applying (18b) for $(X, Y) = (d-1-k, 2)$ implies $P_d(d-2-k, 2) = 0$ for all $k \in \{1, \dots, d-3\}$. Point (e) is therefore straightforward by induction.

For all $d \le k$, the polynomial $Q_{k,d}(X) = P_k(X, d-X)$ satisfies $[Q_{k,d}]_k = [d^k H_k(\frac{2X-d}{d})]_k \ne 0$ for all $k \ge 2$. Therefore, $\deg(Q_{k,d}) = k$ for all $k \ge 0$, and $\{Q_{0,d}, Q_{1,d}, \dots, Q_{d,d}\}$ is a basis of $\mathcal{P}'_d := \{Q \in \mathbb{R}[X] : \deg(Q) \le d\}$. Since $\varphi(Q) = (Q(0), \dots, Q(d))$ defines a linear isomorphism from $\mathcal{P}'_d$ to $\mathbb{R}^{d+1}$, we deduce that $\{\varphi(Q_{0,d}), \dots, \varphi(Q_{d,d})\}$ is a basis of $\mathbb{R}^{d+1}$, which is equivalent to (f).

Point (g) is a simple consequence of points (e) and (f) and of formula (18a) with $X = j$ and $Y = d-j-1$.

Because of point (e) above, (h) is obvious if $k \le d-1$ or $k \le d' - 1$. So let us assume that $d, d' \le k$. Multiplying (18a) by $(X+Y-d+1)$ and applying (18b) to both terms on the l.h.s. yields

$$ (2XY - d(d-1))P_d(X, Y) = XY(P_d(X+1, Y-1) + P_d(X-1, Y+1)). \quad (30) $$

This means that, for all $k \ge 2$ and $2 \le d \le k$, the vector $(P_d(i, k-i))_{1 \le i \le k-1}$ is a right eigenvector of the matrix $A_k = (a_{i,j}^{(k)})_{1 \le i,j \le k-1}$ for the eigenvalue $-d(d-1)$, where

$$ \begin{align*} a_{i,i}^{(k)} &= -2i(k-i) \text{ for } 1 \le i \le k-1, \\ a_{i,i+1}^{(k)} &= i(k-i) \text{ for } 2 \le i \le k-1, \\ a_{i,i-1}^{(k)} &= i(k-i) \text{ for } 1 \le i \le k-2 \\ \text{and } a_{ij}^{(k)} &= 0 \text{ for } |i-j| \ge 2. \end{align*} $$

It is straightforward to check that the matrix $A_k$ is self-adjoint for the inner product $\langle \cdot, \cdot \rangle_\mu$, where $\mu_i = 1/i(k-i)$, which implies that two right eigenvectors of $A_k$ corresponding to different eigenvalues are orthogonal w.r.t. this inner product. This yields (21) for $d \ne d'$.
Finally, fix $2 \le d \le k$. Using (18a) and (18b), we have

$$ 
\begin{align*}
& (k+d) \sum_{i=1}^{k-1} \frac{P_d(i, k-i)^2}{i(k-i)} \\
&= \sum_{i=1}^{k-1} P_d(i, k-i) \frac{P_d(i+1, k-i)}{k-i} + \sum_{i=1}^{k-1} P_d(i, k-i) \frac{P_d(i, k-i+1)}{i} \\
&= \sum_{i=1}^{k} P_d(i, k-i+1) \frac{P_d(i-1, k-i+1)}{k-i+1} + \sum_{i=1}^{k} P_d(i, k-i+1) \frac{P_d(i, k-i)}{i} \\
&= (k-d+1) \sum_{i=1}^{k} \frac{P_d(i, k-i+1)^2}{i(k-i+1)}.
\end{align*}
$$

Applying this equality inductively, we deduce that

$$ \sum_{i=1}^{k-1} \frac{P_d(i, k-i)^2}{i(k-i)} = C \binom{k+d-1}{2d-1} $$

for some constant $C$.

Now, by point (a) and Proposition 2.2, we have

$$ \frac{1}{k^{2d-1}} \sum_{i=1}^{k-1} \frac{P_d^2(i, k-i)}{i(k-i)} \sim \frac{4}{k} \sum_{i=1}^{k-1} \frac{H_d^2\left(\frac{2i-k}{k}\right)}{1 - \left(\frac{2i-k}{k}\right)^2} \to 2 \int_{-1}^{1} \frac{H_d^2(x)}{1-x^2} dx = 2 $$

as $k \to +\infty$. Thus $C = 2$ and the proof of (h) is completed. $\square$

**Proof of Proposition 3.3** In [15], the authors construct a family of functions of two variables satisfying relations close to (18), which they use to study neutral, multitype population processes with non-zero mutation or immigration. These functions are expressed in terms of the Hahn polynomials, defined for fixed parameters $\alpha > -1$, $\beta > -1$ and $N \in \mathbb{N}$ by

$$ Q_d(x; \alpha, \beta, N) = {}_3F_2 \left( \begin{array}{c} -d, -x, d+\alpha+\beta+1 \\ \alpha+1, -N+1 \end{array} ; 1 \right), \qquad (31) $$

for all integer $d \ge 0$. Karlin and McGregor proved that the rational function

$$ \phi_d(X, Y) = Q_d(X; \alpha, \beta, X + Y + 1) = \sum_{k=0}^{d} \frac{(-d)_k (-X)_k (d + \alpha + \beta + 1)_k}{(\alpha + 1)_k (-X - Y)_k k!} $$
satisfies

$$ (X + \alpha + 1)\phi_d(X + 1, Y) + (Y + \beta + 1)\phi_d(X, Y + 1) = (X + Y)\phi_d(X, Y) \quad (32) $$

$$ X\phi_d(X - 1, Y) + Y\phi_d(X, Y - 1) = \frac{(X + Y + 1 - d)(X + Y + d + \alpha + \beta + 2)}{X + Y + 1}\phi_d(X, Y). \quad (33) $$

Let us define

$$ \begin{aligned} \psi_d(X, Y) &= (-X - Y)_d \lim_{\beta \to -1} \lim_{\alpha \to -1} (\alpha + 1)\phi_d(X, Y) \\ &= (-X - Y)_d \sum_{k=1}^{d} \frac{(-d)_k (-X)_k (d-1)_k}{(k-1)!(-X-Y)_k!} \end{aligned} $$

Passing to the limit in (32) and (33) proves that $\psi_d$ satisfies (18). Since $\psi_d$ is a polynomial, Theorem 3.1 entails (22). It only remains to check that $C_d\psi_d(X, Y)$ has its $(d-1, 1)$ coefficient equal to $(-1)^{d+2}\sqrt{d(d-1)(2d-1)}$. The $(d-1, 1)$ coefficient of $\psi_d(X, Y)$ is

$$ \begin{aligned} \sum_{k=1}^{d} \frac{(-d)_k (d-1)_k}{(k-1)!k!} (d-k) &= (-1)^{d+1} d(d-1)^2 \sum_{k=1}^{d} \frac{(-d+2)_{k-1}(d)_{k-1}}{(k-1)!k!} \\ &= (-1)^{d+1} d(d-1) {}_2F_1\left(-\begin{array}{c} d+2, d \\ 2 \end{array}; 1\right). \end{aligned} $$

The Chu-Vandermonde formula (cf. e.g. [9]) implies that

$$ {}_2F_1\left(-\begin{array}{c} d+2, d \\ 2 \end{array}; 1\right) = (-1)^d \frac{(d-2)!}{(d-1)!} = \frac{(-1)^d}{d-1}, $$

which gives the expression of $C_d$.

<figure><img src="image_2.png" /></figure>

# 4 Spectral decomposition of neutral two-dimensional Markov chains

In this section, we consider neutral extensions of the two-dimensional birth and death chains in $\mathbb{Z}^2$ described in the introduction. In the sequel, this family will be called N2dMC, for neutral two-dimensional Markov chains.
A N2dMC $(X_t, Y_t)_{t \in \mathbb{Z}_+}$ is constructed by specifying first the Markov dynamics of $Z_t = X_t + Y_t$ in $\mathbb{Z}_+$. Assume that its transition matrix is

$$ \Pi_0 = (p_{n,m})_{n,m \ge 0}, $$

where $\sum_{m \ge 0} p_{n,m} = 1$ for all $n \ge 1$ and the state 0 is absorbing ($p_{0,0} = 1$). Then, the process $(X_t, Y_t)_{t \in \mathbb{Z}_+}$ is constructed as follows: if there is a birth at time $t$ (i.e. if $Z_{t+1} > Z_t$), the types of the new individuals are *successively* picked at random in the population; if there is a death at time $t$ (i.e. if $Z_{t+1} < Z_t$), the types of the killed individuals are *successively* picked at random in the population; finally, if $Z_{t+1} = Z_t$, then $X_{t+1} = X_t$ and $Y_{t+1} = Y_t$.

For example, the transition probability from $(i, j)$ to $(i + k, j + l)$ for $(i, j) \in \mathbb{Z}_+^2$ and $k, l \ge 0$, $k + l \ge 1$ is

$$ \binom{k+l}{k} \frac{i(i+1)\dots(i+k-1)j(j+1)\dots(j+l-1)}{(i+j)(i+j+1)\dots(i+j+k+l-1)} p_{i+j,i+j+k+l}. $$

After some algebra, one gets the following formulas for the transition probabilities: for all $l \ge 0$ and $k \ge 0$ such that $l + k \ge 1$, the Markov chain $(X_n, Y_n)_{n \ge 0}$ has transitions from $(i, j) \in \mathbb{Z}_+^2$ to

$$ \begin{align} (i+k, j+l) \text{ w. p. } \pi_{(i,j),(i+k,j+l)} &:= \frac{\binom{i+k-1}{k} \binom{j+l-1}{l}}{\binom{i+j+k+l-1}{k+l}} p_{i+j, i+j+k+l} \tag{34} \\ (i-k, j-l) \text{ w. p. } \pi_{(i,j),(i-k,j-l)} &:= \frac{\binom{i}{k} \binom{j}{l}}{\binom{i+j}{k+l}} p_{i+j, i+j-k-l} \notag \\ (i,j) \text{ w. p. } \pi_{(i,j),(i,j)} &:= p_{i+j, i+j}, \notag \end{align} $$

with the convention that $\binom{i}{j} = 0$ if $i < 0$, $j < 0$ or $j > i$. In particular, once one component of the process is 0, it stays zero forever. We denote by

$$ \Pi := (\pi_{(i,j),(k,l)})_{(i,j),(k,l) \in \mathbb{Z}_+^2} $$

the transition matrix of the Markov process $(X, Y)$.

The state space of the Markov chain $Z$ will be denoted by $S_Z$, and the state space of $(X, Y)$ by $S$. We are going to consider two cases: the case where $Z$ has finite state space $S_Z := \{0, 1, \dots, N\}$ for some $N \ge 0$, and the case where $Z$ has infinite state space $S_Z := \mathbb{Z}_+$. In the first case, the state space of $(X, Y)$ is the set

$$ \mathcal{T}_N := \{(i, j) \in \mathbb{Z}_+^2 : i + j \le N\}. \tag{35} $$
In the second case, $S := \mathbb{Z}_+^2$.

We also define the sets $S^* := S \cap \mathbb{N}^2$, $\mathcal{T}_N^* := \mathcal{T}_N \cap \mathbb{N}^2$ and $S_Z^* := S_Z \setminus \{0\}$. Finally, let $\tilde{\Pi}_0$ be the restriction of the matrix $\Pi_0$ to $S_Z^*$ (i.e. the matrix obtained from $\Pi_0$ by suppressing the row and column of index 0) and let $\tilde{\Pi}$ be the restriction of the matrix $\Pi$ to $S^*$.

Extending the usual definition for Markov chains, we say that a matrix $M = (m_{ij})_{i,j \in \mathcal{A}}$ is reversible with respect to the measure $(\mu_i)_{i \in \mathcal{A}}$ if $\mu_i > 0$ and $\mu_i m_{ij} = \mu_j m_{ji}$ for all $i, j \in \mathcal{A}$.

For all $d \ge 0$, we define

$$ \mathcal{V}_d := \{ v \in \mathbb{R}^S : v_{i,j} = P_d(i,j) u_{i+j} \text{ with } u \in \mathbb{R}^{S_Z} \}, $$

where we recall that the polynomials $P_d$ are defined in Theorem 3.1 and $P_1 = P_1^{(1)}$. Note that, by Proposition 3.2 (e), a vector $v_{i,j} = P_d(i,j)u_{i+j} \in \mathcal{V}_d$ is characterized by the values of $u_k$ for $k \ge d$ only.

For all $d \ge 0$, we also define the matrix $\Pi_d := (p_{n,m}^{(d)})_{(n,m) \in S, n \ge d, m \ge d}$, where for all $(n, m) \in S$ such that $n \ge d$ and $m \ge d$,

$$ p_{n,m}^{(d)} := \begin{cases} \frac{\binom{m+d-1}{m-n}}{\binom{m-1}{m-n}} p_{n,m} & \text{if } m > n, \\ \frac{\binom{n-d}{n-m}}{\binom{n}{n-m}} p_{n,m} & \text{if } m < n, \\ p_{n,n} & \text{if } m = n. \end{cases} \qquad (36) $$

All these notation, as well as those introduced in the rest of the paper, are gathered for convenience in Appendix A.

The following result is the basis of all results in this section.

**Proposition 4.1** For all $d \ge 0$, the vector space $\mathcal{V}_d$ is stable for the matrix $\Pi$. In addition, for all $v_{i,j} = P_d(i,j)u_{i+j} \in \mathcal{V}_d$,

$$ (\Pi v)_{i,j} = P_d(i,j) (\Pi_d u)_{i+j}. $$

**Proof** Using (18a) inductively, we have

$$ \sum_{k=0}^{n} \binom{n}{k} X(X+1) \dots (X+k-1)Y(Y+1) \dots (Y+n-k-1)P_d(X+k, Y+n-k) \\ = (X+Y+d)(X+Y+d+1) \dots (X+Y+d+n-1)P_d(X,Y) $$
for all $d \ge 0$ and $n \ge 1$, which can be written as

$$ \sum_{k=0}^{n} \binom{X+k-1}{k} \binom{Y+n-k-1}{n-k} P_d(X+k, Y+n-k) \\ = \binom{X+Y+d+n-1}{n} P_d(X,Y) $$

for $(X, Y) \in \mathbb{Z}_+^2$. Similarly, an inductive use of (18b) yields

$$ \sum_{k=0}^{n} \binom{X}{k} \binom{Y}{n-k} P_d(X-k, Y-n+k) = \binom{X+Y-d}{n} P_d(X,Y) $$

for all $d \ge 0$, $n \ge 1$ and $(X, Y) \in \mathbb{Z}_+^2$, where we recall the convention $\binom{a}{b} = 0$ if $a < 0$, $b < 0$ or $b > a$.

Proposition 4.1 then easily follows from these equations and from the transition probabilities (34). $\square$

# 4.1 The case of finite state space

## 4.1.1 Eigenvectors of $\Pi$ for finite state spaces

In the case where $Z$ has finite state space, the main result of this section is the following.

**Theorem 4.2** Assume that $S_Z = \mathcal{T}_N$ for some $N \ge 0$.

(a) For all $d \ge 0$ and all right eigenvector $(u_n)_{n \in S_Z, n \ge d}$ of $\Pi_d$ for some eigenvalue $\theta$, the vector

$$ v_{(i,j)} = P_d(i,j)u_{i+j}, \quad (i,j) \in S \qquad (37) $$

is a right eigenvector of the matrix $\Pi$ for the same eigenvalue, where the polynomials $P_d$ are defined in Theorem 3.1 and where $P_1 = P_1^{(1)}$. In addition, if $d \ge 2$, $v$ is also a right eigenvector of the matrix $\tilde{\Pi}$ for the same eigenvalue.

(b) All the right eigenvectors of $\Pi$ are of the form (37), or possibly a linear combination of such eigenvectors in the case of multiple eigenvalues.
(c) Assume that $\tilde{\Pi}_0$ admits a positive reversible measure $(\mu_n)_{n \in S_Z^*}$. Then, the matrix $\tilde{\Pi}$ is reversible w.r.t. the measure

$$ \nu_{(i,j)} := \frac{(i+j)\mu_{i+j}}{ij}, \quad \forall (i,j) \in S^*, \qquad (38) $$

and hence is diagonalizable in a basis orthonormal for this measure, composed of vectors of the form (37) for $d \ge 2$.

In addition, $\Pi$ is diagonalizable in a basis of eigenvectors of the form (37) for $d \ge 0$.

Hence, the right eigenvectors of the transition matrix of a finite N2dMC can be decomposed as the product of two terms, one depending on each population size, but “universal” in the sense that it does not depend on the transitions matrix $\Pi_0$ of $Z$, and the other depending on the matrix $\Pi_0$, but depending only on the total population size.

**Remark 4.3** There is some redundancy among the right eigenvectors of $\Pi$ of the form $P(i, j)u_{i+j}$ for $P = P_0$, $P_1^{(1)}$ or $P_1^{(2)}$: if $u$ is a right eigenvector of $\Pi_1$, the vectors

$$ (iu_{i+j})_{(i,j) \in S} \quad \text{and} \quad (ju_{i+j})_{(i,j) \in S} $$

are right eigenvectors of $\Pi$ for the same eigenvalue. In particular, $iu_{i+j} + ju_{i+j}$ is an eigenvector of $\Pi$ of the form $P_0(i, j)u'_{i+j}$. This will also be true when $Z$ has infinite state space.

**Remark 4.4** In the following proof (and also in the case of infinite state space), no specific use is made of the fact that the matrix $\Pi$ is stochastic. Therefore, Theorem 4.2 also holds true in the case of a continuous-time N2dMC, where the matrix $\Pi_0$ is now the infinitesimal generator of the process $Z_t = X_t + Y_t$.

**Proof** Point (a) is an easy consequence of Proposition 4.1.

For all $0 \le d \le N$, the matrix $\Pi_d$ is conjugate with its Jordan normal form. Let $\{u^{(d),k}\}_{d \le k \le N}$ denote the basis of $\mathbb{C}^{N-d+1}$ corresponding to this normal form, where $u^{(d),k} = (u_d^{(d),k}, u_{d+1}^{(d),k}, \dots, u_N^{(d),k})$. Then, the family of vectors

$$ \mathcal{F} := \{ (P_d(i, j)u_{i+j}^{(d),k})_{(i,j) \in \mathcal{T}_N} : 0 \le d \le N, d \le k \le N \} $$
is composed of $(N+1)+N+(N-1)+...+1 = (N+1)(N+2)/2 = |\mathcal{T}_N|$ elements. Moreover, one can prove that it is linearly independent as follows: since $\{u^{(d),k}\}_{d \le k \le N}$ is a basis of $\mathbb{C}^{N-d+1}$, it is sufficient to check that

$$ \sum_{d=0}^{N} P_d(i, j) v_{i+j}^{(d)} = 0, \quad \forall (i, j) \in \mathcal{T}_N \qquad (39) $$

implies that $v^{(d)} = 0$ for all $0 \le d \le N$, where $v^{(d)} = (v_d^{(d)},...,v_N^{(d)}) \in \mathbb{C}^{N-d+1}$, seen as a subspace of $\mathbb{C}^{N+1}$ by putting the $d$ first coordinates to be zero. Given $k \le N$, the equality (39) for $i + j = k$ combined with Proposition 3.2 (f) yields $v_k^{(0)} = ... = v_k^{(k)} = 0$.

Therefore, $\mathcal{F}$ is a basis of $\mathbb{C}^{\mathcal{T}_N}$ and, by point (a), the matrix $\Pi$ has a Jordan normal form in this basis. Point (b) is then straightforward.

If $\tilde{\Pi}_0$ admits a positive reversible measure $\mu$, it is straightforward to check that the vector $\nu$ in (38) is a reversible measure for $\tilde{\Pi}$, and hence the first part of Point (c) is true.

In addition, the matrix $\Pi_1$ is reversible w.r.t. the measure

$$ \mu_n^{(1)} := 2 n^2 \mu_n, \quad n \in S_Z^*, $$

which implies that $\Pi_1$ admits a basis of right eigenvectors orthonormal w.r.t. $\mu^{(1)}$. Similarly, $\tilde{\Pi}_0$ admits a basis of right eigenvectors orthonormal w.r.t. $\mu$. By Point (a), this gives $N+1$ (resp. $N$) right eigenvectors of $\Pi$ of the form (37) for $d = 0$ (resp. $d = 1$). Together with the basis of right eigenvectors of $\tilde{\Pi}$ obtained above (extended by zero on $\{0\} \times \mathbb{N}$ and $\mathbb{N} \times \{0\}$), this gives a basis of eigenvectors of $\Pi$ and ends the proof of (c). $\square$

### 4.1.2 Example: 3-colors urn model (or 3-types Moran model)

The class of transition matrices given in (34) can be obtained by composition and linear combinations of the transition matrices $\Pi^{(n)+} = (\pi_{(i,j),(k,l)}^{(n)+})_{(i,j),(k,l) \in S}$, $\Pi^{(n)-} = (\pi_{(i,j),(k,l)}^{(n)-})_{(i,j),(k,l) \in S}$ and $\Pi^{(n)} = (\pi_{(i,j),(k,l)}^{(n)})_{(i,j),(k,l) \in S}$, where for all
$n \ge 1$

$$ \pi^{(n)+} = \begin{cases} \frac{i}{i+j}, & \text{if } k = i+1, l = j, i+j = n \\ \frac{j}{i+j}, & \text{if } k = i, l = j+1, i+j = n \\ 0, & \text{otherwise,} \end{cases} $$

$$ \pi^{(n)-} = \begin{cases} \frac{i}{i+j}, & \text{if } k = i-1, l = j, i+j = n \\ \frac{j}{i+j}, & \text{if } k = i, l = j-1, i+j = n \\ 0, & \text{otherwise,} \end{cases} $$

and for all $n \ge 0$

$$ \pi^{(n)} = \begin{cases} 1, & \text{if } k = i, l = j, i+j = n \\ 0, & \text{otherwise.} \end{cases} $$

One easily checks, first that the vector spaces $V_d$ for all $d \ge 0$ are stable for all these matrices, and second that, for the matrix (34),

$$ \Pi = \sum_{n<m \in S_Z} p_{n,m} (\Pi^{(n)+})^{m-n} + \sum_{n>m \in S_Z} p_{n,m} (\Pi^{(n)-})^{n-m} \sum_{n \in S_Z} p_{n,m} \Pi^{(n)}. $$

Hence, the vector spaces $V_d$ are trivially stable for such matrices.

One may however recover a much larger class of matrices for which $V_d$ are stable vector spaces, by considering the algebra of matrices spanned by the matrices $\Pi^{(n)\pm}$ and $\Pi^{(n)}$. Below, we study in detail such an example.

Consider an urn with $N$ balls of three different colors and consider the following process: one picks a ball at random in the urn, notes its color, puts it back in the urn, picks another ball in the urn and replaces it by a ball of the same color as the first one. The number of balls of each colors then forms a Markov chain, which can be viewed as the embedded Markov chain of the 3-types Moran model, defined as follows: consider a population of $N$ individuals, with 3 different types. For each pair of individuals, at rate 1, the second individual is replaced by an individual of the same type as the first individual in the pair.

Let $i$ denote the number of balls of the first color, and $j$ the number of balls of the second color. Then, there are $N - i - j$ balls of the third color. The transition probabilities of this Markov chain are as follows: given that
the current state of the process is $(i, j)$, the state at the next time step is

$$
\begin{align*}
(i+1, j) \quad & \text{with probability} \quad \frac{i(N-i-j)}{N^2}, \\
(i-1, j) \quad & \text{with probability} \quad \frac{i(N-i-j)}{N^2}, \\
(i, j+1) \quad & \text{with probability} \quad \frac{j(N-i-j)}{N^2}, \\
(i, j-1) \quad & \text{with probability} \quad \frac{j(N-i-j)}{N^2}, \\
(i+1, j-1) \quad & \text{with probability} \quad \frac{ij}{N^2}, \\
(i-1, j+1) \quad & \text{with probability} \quad \frac{ij}{N^2}, \\
(i, j) \quad & \text{with probability} \quad \frac{i^2+j^2+(N-i-j)^2}{N^2}.
\end{align*}
$$

These transition probabilities do not not have the form (34). However, a variant of Proposition 4.1, and hence of Theorem 4.2, apply to this process, because of the following observation: let us construct the matrices $\Pi^+$, $\Pi^-$, $\hat{\Pi}^+$, $\tilde{\Pi}^+$ from the matrices $\Pi_0^+ = (p_{n,m}^+ )_{n,m \in S_Z}$, $\Pi_0^- = (p_{n,m}^- )_{n,m \in S_Z}$, $\hat{\Pi}_0^+ = (\hat{p}_{n,m}^+ )_{n,m \in S_Z}$, $\tilde{\Pi}_0^+ = (\tilde{p}_{n,m}^+ )_{n,m \in S_Z}$ respectively, exactly as $\Pi$ was constructed from $\Pi_0$ in (34), where

$$
\begin{align*}
p_{n,n+1}^{+} &= \frac{k}{N}, & p_{n,n}^{+} &= 1 - \frac{k}{N}, & p_{n,m}^{+} &= 0 \text{ otherwise,} \\
p_{n,n-1}^{-} &= \frac{k}{N}, & p_{n,n}^{+} &= 1 - \frac{k}{N}, & p_{n,m}^{+} &= 0 \text{ otherwise,} \\
\hat{p}_{n,n+1}^{+} &= \frac{k(N-k)}{N^2}, & \hat{p}_{n,n}^{+} &= 1 - \frac{k(N-k)}{N^2}, & \hat{p}_{n,m}^{+} &= 0 \text{ otherwise,} \\
\hat{p}_{n,n+1}^{-} &= \frac{k(N-k-1)}{N^2}, & \hat{p}_{n,n}^{+} &= 1 - \frac{k(N-k-1)}{N^2}, & \hat{p}_{n,m}^{+} &= 0 \text{ otherwise.}
\end{align*}
$$

Then the transition matrix of the 3-colors urn model is given by

$$
\Pi = \Pi^{+} \Pi^{-} + \hat{\Pi}^{+} - \tilde{\Pi}^{-}.
$$

In particular, the vector spaces $V_d$ for $0 \le d \le N$ are all stable for this
matrix.

The transition matrix Π has absorbing sets { (i, 0) : 0 ≤ i ≤ N }, { (0, i ) :
0 ≤ i ≤ N } and { ( i, N − i ) : 0 ≤ i ≤ N }, and absorbing states (0, 0), ( N, 0)
and (0, N ). The restriction of the matrix Π on the set

$$
\mathcal{S}^{**} := \{(i,j) : i \ge 1, j \ge 1, i + j \le N - 1\}
$$

admits the reversible measure

$$
\nu_{(i,j)} = \frac{1}{ij(N - i - j)}.
$$

Hence the matrix Π admits a family of right eigenvectors null on the absorb-
ing sets, which forms an orthonormal basis of L²(S**, ν).
One easily checks, using (30), that $v_{(i,j)} = P_d(i,j)u_{i+j}$ is a right eigenvector of $\Pi$ for an eigenvalue $\theta$ if and only if for all $d \le k \le N$

$$ \begin{aligned} \theta'u_k &= (N-k)[(k+d)u_{k+1} - 2ku_k + (k-d)u_{k-1}] \\ &= (N-k)[k(u_{k+1} - 2ku_k + u_{k-1}) + d(u_{k+1} - u_{k-1})], \end{aligned} \quad (40) $$

where

$$ \theta = 1 + \frac{\theta' - d(d-1)}{N^2}. $$

Now, the Hahn polynomials $Q_n(x; \alpha, \beta, N)$ introduced in (31) satisfy (cf. e.g. [15])

$$ \begin{aligned} -n(n+\alpha+\beta+1)Q_d(x) = x(N+\beta-x)Q_n(x-1) &+ (N-1-x)(\alpha+1+x)Q_n(x+1) \\ & - [x(N+\beta-x) + (N-1-x)(\alpha+1+x)]Q_n(x). \end{aligned} $$

Hence, for all $0 \le n \le N-d$, (40) admits the polynomial (in $k$) solution of degree $n$

$$ u_k = Q_n(k-d; 2d-1, -1, N-d+1). $$

If $n \ge 1$, this polynomial must be divisible by $(N-k)$, so we can define the polynomial $R_n^{(N,d)}(X)$ of degree $n-1$ as

$$ (N-X)R_n^{(N,k)}(X) = Q_n(X-d; 2d-1, -1, N-d+1). $$

Obviously, the family of vectors $(1, \dots, 1)$ and $((N-d)R_n^{(N,d)}(d), \dots, R_n^{(N,d)}(N-1), 0)$ for $1 \le n \le N-d$ is linearly independent and hence forms a basis of the vector space $\mathbb{R}^{N-d+1}$ of real vectors indexed by $d, d+1, \dots, N$. In addition, (40) cannot admit any other linearly independent solution and hence, necessarily, $R_n^{(N,d)}(k) = 0$ for all $n > N-d$ and $d \le k \le N-1$.

We have obtained a basis of right eigenvectors of $\Pi$ of the form

$$ \{P_d(i, j)\}_{0 \le d \le N} \cup \{P_d(i, j)(N - i - j)R_n^{(N,d)}(i + j)\}_{0 \le d \le N-1, 1 \le n \le N-d}, $$

and the eigenvalue corresponding to the eigenvector $P_d(i, j)$ if $n = 0$, or $P_d(i, j)(N - i - j)R_n^{(N,d)}(i + j)$ if $n \ge 1$, is

$$ \theta_{d,n} := 1 - \frac{n(n-1) + 2nd - d(d-1)}{N^2} = 1 - \frac{(d+n)(d+n-1)}{N^2}. $$
Similarly as in the proof of Proposition 3.3, this family of eigenvectors can be seen as a singular limit case of those obtained in [15] for the multitype Moran model with mutation or immigration.

Note that in the case of the 2-colors urn model, one can easily check that a basis of right eigenvectors of the corresponding transition matrix is given by $(1, \dots, 1)$ and $(NR_n^{(N,0)}(0), \dots, R_n^{(N,0)}(N-1), 0)$ for $1 \le n \le N$. Hence the spectrum is the same in the 2- and 3-colors urn models, although the multiplicity of each eigenvalue is different. In the case of two colors, the eigenvalues have the form $1 - k(k-1)/N^2$ for $0 \le k \le N$, each with multiplicity 1 (except for the eigenvalue 1, with multiplicity 2). In the case of three colors, the eigenvalue $1 - k(k-1)/N^2$ has multiplicity $k+1$ (except for the eigenvalue 1, which has multiplicity 3).

Concerning the eigenvectors in $\mathbb{L}^2(S^{**}, \nu)$, they are given by

$$ \{ij(N - i - j)Q_d(i, j)R_n^{(N,d)}(i + j)\}_{2 \le d \le N-1, 1 \le n \le N-d}, $$

and for all $3 \le k \le N$, the eigenspace for the eigenvalue $1 - k(k-1)/N^2$ is

$$ V_k := \mathrm{Vect}\{ij(N - i - j)Q_2(i, j)R_{k-2}^{(N,2)}(i + j), \dots, \\ ij(N - i - j)Q_{k-1}(i, j)R_1^{(N,k-1)}(i + j)\}. $$

We shall end the study of this example by giving an apparently non-trivial relation between the polynomials $P_d$ and $R_n^{(N,d)}$. Because of the symmetry of the colors, we have

$$ \begin{aligned} V_k &= \mathrm{Vect}\{ij(N - i - j)Q_2(i, N - i - j)R_{k-2}^{(N,2)}(N - j), \dots, \\ & \qquad ij(N - i - j)Q_{k-1}(i, N - i - j)R_1^{(N,k-1)}(N - j)\} \\ &= \mathrm{Vect}\{ij(N - i - j)Q_2(N - i - j, j)R_{k-2}^{(N,2)}(N - i), \dots, \\ & \qquad ij(N - i - j)Q_{k-1}(N - i - j, j)R_1^{(N,k-1)}(N - i)\}, \end{aligned} $$

and hence

$$ \begin{aligned} \mathrm{Vect}\{Q_2(i, j)R_{n-2}^{(N,2)}(i + j), \dots, Q_{n-1}(i, j)R_1^{(N,n-1)}(i + j)\}_{(i,j) \in S^{**}} \\ = \mathrm{Vect}\{Q_2(i, N - i - j)R_{n-2}^{(N,2)}(N - j), \dots, Q_{n-1}(i, N - i - j)R_1^{(N,n-1)}(N - j)\}_{(i,j) \in S^{**}} \\ = \mathrm{Vect}\{Q_2(N - i - j, j)R_{n-2}^{(N,2)}(N - i), \dots, Q_{n-1}(N - i - j, j)R_1^{(N,n-1)}(N - i)\}_{(i,j) \in S^{**}}. \end{aligned} $$
## 4.2 The case of infinite state space

The goal of this subsection is to extend Theorem 4.2 to the case where $Z$ has infinite state space and the matrix $\Pi_0$ is a compact operator and admits a reversible measure. To this aim, we need first some approximation properties of $\Pi$ by finite rank operators.

### 4.2.1 Approximation properties

Recall that $\Pi_0$ is a Markov kernel on $\mathbb{Z}_+$ absorbed at 0 and that $\tilde{\Pi}_0$ denotes its restriction to $\mathbb{N}$ (equivalently, $\tilde{\Pi}_0$ is the sub-Markovian kernel on $\mathbb{N}$ corresponding to $\Pi_0$ with a Dirichlet condition at 0). We assume that $\tilde{\Pi}_0$ is reversible with respect to a positive measure $\mu$ on $\mathbb{N}$ (not necessarily finite). For any $N \in \mathbb{N}$, consider the sub-Markovian kernel $\tilde{\Pi}_0^{(N)}$ on $\mathbb{N}$ defined by

$$ \forall x, y \in \mathbb{N}, \quad \tilde{\Pi}_0^{(N)}(x, y) := \begin{cases} \Pi_0(x, y) & \text{if } x, y \le N \\ 0 & \text{otherwise.} \end{cases} $$

In other words, $\tilde{\Pi}_0^{(N)} = \text{Pr}_N \tilde{\Pi}_0 \text{Pr}_N$, where $\text{Pr}_N$ is the projection operator defined by $\text{Pr}_N(u_1, u_2, \dots) = (u_1, \dots, u_N, 0, \dots)$.

The kernel $\tilde{\Pi}_0^{(N)}$ is not Markovian (since $\tilde{\Pi}_0^{(N)}(x, \mathbb{N}) = 0$ for $x > N$), but it can be seen as the restriction to $\mathbb{N}$ of a unique Markovian kernel $\Pi_0^{(N)}$ on $\mathbb{Z}_+$ absorbed at 0. With this interpretation, we can construct the matrix $\tilde{\Pi}^{(N)}$ from $\Pi_0^{(N)}$ exactly as the matrix $\tilde{\Pi}$ was constructed from $\Pi_0$ in the beginning of Section 4.

Of course $\tilde{\Pi}_0^{(N)}$ remains reversible with respect to $\mu$, and thus, like $\tilde{\Pi}_0$, it can be extended into a self-adjoint operator on $\mathbb{L}^2(\mathbb{N}, \mu)$. We denote by $\|\cdot\|_0$ the natural operator norm on the set of bounded operators on $\mathbb{L}^2(\mathbb{N}, \mu)$, namely if $K$ is such an operator,

$$ \|K\|_0 := \sup_{u \in \mathbb{L}^2(\mathbb{N}, \mu) \setminus \{0\}} \frac{\|Ku\|_{\mu}}{\|u\|_{\mu}}. $$

If furthermore $K$ is self-adjoint, we have, via spectral calculus,

$$ \|K\|_0 := \sup_{u \in \mathbb{L}^2(\mathbb{N}, \mu) \setminus \{0\}} \frac{|\langle u, Ku \rangle_{\mu}|}{\|u\|_{\mu}^2}. $$

The next result gives a simple compactness criterion for $\tilde{\Pi}_0$.
**Lemma 4.5** The operator $\tilde{\Pi}_0$ acting on $L^2(\mathbb{N}, \mu)$ is compact if and only if

$$ \lim_{N \to \infty} \| \tilde{\Pi}_0 - \tilde{\Pi}_0^{(N)} \|_0 = 0. $$

**Proof** Since for any $N \in \mathbb{Z}_+$, $\tilde{\Pi}_0^{(N)}$ has finite range, if the above convergence holds, $\tilde{\Pi}_0$ can be strongly approximated by finite range operators and it is well-known that this implies that $\tilde{\Pi}_0$ is compact.

The converse implication can be proved adapting a standard argument for compact operators: assume that $\tilde{\Pi}_0$ is compact and let $\varepsilon > 0$ be fixed. Then $\tilde{\Pi}_0(B)$ is compact, where $B$ is the closed unit ball of $L^2(\mathbb{N}, \mu)$. Hence

$$ \tilde{\Pi}_0(B) = \bigcup_{i=1}^{n} B(\psi^{(i)}, \varepsilon), $$

for some $n < +\infty$ and $\psi^{(1)}, \dots, \psi^{(n)} \in \tilde{\Pi}_0(B)$, where $B(\psi, \varepsilon)$ is the closed ball centered at $\psi$ with radius $\varepsilon$. For any $i \le n$, since $\psi^{(i)} \in L^2(\mathbb{N}, \mu)$, there exists $N_i$ such that $\sum_{k \ge N_i} (\psi_k^{(i)})^2 \mu_k \le \varepsilon$, and thus

$$ \tilde{\Pi}_0(B) = \bigcup_{i=1}^{n} B(\mathrm{Pr}_{N_i} \psi^{(i)}, 2\varepsilon). $$

In other words, for all $\varphi \in B$, there exists $i \le n$ such that $\|\tilde{\Pi}_0\varphi - \mathrm{Pr}_{N_i}\psi^{(i)}\|_{\mu} \le 2\varepsilon$. This implies that

$$ \|\tilde{\Pi}_0\varphi - \mathrm{Pr}_N \tilde{\Pi}_0\varphi\|_{\mu} \le 2\varepsilon, $$

where $N = \sup\{N_1, \dots, N_n\}$, i.e. $\|\tilde{\Pi}_0 - \mathrm{Pr}_N \tilde{\Pi}_0^{(N)}\|_0 \le 2\varepsilon$. Since $\tilde{\Pi}_0^{(N)} = \mathrm{Pr}_N \tilde{\Pi}_0 \mathrm{Pr}$, we obtain that

$$ \lim_{N \to +\infty} \| \tilde{\Pi}_0 \mathrm{Pr}_N - \tilde{\Pi}_0^{(N)} \|_0 = 0. $$

In order to complete the proof, it only remains to check that

$$ \lim_{N \to +\infty} \| \tilde{\Pi}_0 - \tilde{\Pi}_0 \mathrm{Pr}_N \|_0 = 0. $$

If this was false, one could find a sequence $(\varphi^{(N)})_{N \ge 1}$ in $B$ such that $\varphi_k^{(N)} = 0$ for all $k \le N$ and $\|\tilde{\Pi}_0\varphi^{(N)}\|_{\mu}$ would not converge to 0. Such a sequence $(\varphi^{(N)})_{N \ge 1}$ weakly converges to 0. Now, another usual characterization of
compact operators is the fact that the image of weakly converging sub-sequences strongly converges to the image of the limit. In other words, $\|\tilde{\Pi}_0\varphi^{(N)}\|_{\mu} \to 0$. This contradiction ends the proof of Lemma 4.5. $\square$

The interest of $\tilde{\Pi}_0^{(N)}$ is that it brings us back to the finite situation. Let $\hat{\Pi}_0^{(N)}$ be the restriction of $\tilde{\Pi}_0^{(N)}$ to $[1, N]$, which can be seen as a $N \times N$ matrix. We have for instance that the spectrum of $\tilde{\Pi}_0^{(N)}$ is the spectrum of $\hat{\Pi}_0^{(N)}$ plus the eigenvalue 0.

We are now going to see how the results of Section 4.1 are affected by the change from $\Pi_0$ to $\Pi_0^{(N)}$. More generally, we consider two Markov kernels $\Pi_0$ and $\Pi_0'$ on $\mathbb{Z}_+$ absorbed at 0, whose restrictions to $\mathbb{N}$, $\tilde{\Pi}_0$ and $\tilde{\Pi}_0'$, are both reversible with respect to $\mu$. We associate to them $\Pi$ and $\Pi'$ defined on $\mathbb{Z}_+^2$ as in (34), and their respective restriction to $\mathbb{N}^2$, $\tilde{\Pi}$ and $\tilde{\Pi}'$. We also define the matrices $\Pi_d$ and $\Pi_d'$ for $d \ge 1$, as in (36). Note that $\tilde{\Pi}$ and $\tilde{\Pi}'$ are reversible with respect to $\nu$, defined in (38) and it is straightforward to check that, for any $d \ge 1$, $\Pi_d$ and $\Pi_d'$ are both reversible w.r.t. $\mu^{(d)} = (\mu_n^{(d)})_{n \in S_d, n \ge d}$ defined by

$$ \mu_n^{(d)} := 2n \binom{n+d-1}{2d-1} \mu_n, \quad n \in \mathbb{N}, \quad n \ge d. \qquad (41) $$

We will denote $\|\cdot\|$ and $\|\cdot\|_d$ the operator norms in $\mathbb{L}^2(\mathbb{N}^2, \nu)$ and $\mathbb{L}^2(\mathbb{N}_d, \mu^{(d)})$, where $\mathbb{N}_d := \{d, d+1, \dots\}$. The next result shows that, if one takes $\Pi_0' = \Pi_0^{(N)}$, the approximation of $\tilde{\Pi}_0$ by $\tilde{\Pi}_0^{(N)}$ behaves nicely.

**Proposition 4.6** We always have

$$ \|\tilde{\Pi} - \tilde{\Pi}'\| = \sup_{d \ge 2} \|\Pi_d - \Pi_d'\|_d $$

Furthermore, if $\tilde{\Pi}_0 - \tilde{\Pi}_0' \ge 0$ (in the sense that all the entries of this infinite matrix are non-negative), then

$$ \forall d \ge 1, \quad \|\Pi_{d+1} - \Pi_{d+1}'\|_{d+1} \le \|\Pi_d - \Pi_d'\|_d \qquad (42) $$

and

$$ \|\Pi_1 - \Pi_1'\|_1 = \|\tilde{\Pi}_0 - \tilde{\Pi}_0'\|_0. \qquad (43) $$

In particular,

$$ \|\tilde{\Pi} - \tilde{\Pi}'\| = \|\Pi_2 - \Pi_2'\|_2 \le \|\tilde{\Pi}_0 - \tilde{\Pi}_0'\|_0. $$
**Proof** For $d \ge 2$, denote by $V'_d$ the set of $v \in L^2(\mathbb{N}^2, \nu)$ of the form

$$ \forall i, j \in \mathbb{N}, \quad v_{i,j} = P_d(i, j)u_{i+j} $$

with $u \in L^2(\mathbb{N}_d, \mu^{(d)})$. We denote by $v(d, u)$ the sequence $v$ defined by the above r.h.s. The definitions of $\nu$, $\mu^{(d)}$ and Proposition 3.2 (h) enable us to see that the mapping

$$ L^2(\mathbb{N}_d, \mu^{(d)}) \ni u \mapsto v(d, u) \in L^2(\mathbb{N}^2, \nu) $$

is an isometry.

Proposition 3.2 (h) also shows that $V'_d$ and $V'_d$ are orthogonal subspaces of $L^2(\mathbb{N}^2, \nu)$ for all $d, d' \ge 2$, $d \ne d'$. We actually have

$$ L^2(\mathbb{N}^2, \nu) = \bigoplus_{d \ge 2} V'_d. \qquad (44) $$

Indeed, let $v \in L^2(\mathbb{N}^2, \nu)$ be orthogonal to $V'_d$ for all $d \ge 2$. For all $d \ge 2$, we define the vector

$$ v_l^{(d)} = \frac{1}{\binom{l+d-1}{2d-1}} \sum_{i=1}^{l-1} v_{i,l-i} \frac{P_d(i, l-i)}{i(l-i)}, \quad l \ge d. $$

The Cauchy-Schwartz inequality and Proposition 3.2 (h) imply that $v^{(d)} \in L^2(\mathbb{N}_d, \mu^{(d)})$ and since $v$ is orthogonal to $V'_d$, by the definition of $V'_d$, the vector $v^{(d)}$ is orthogonal to $L^2(\mathbb{N}_d, \mu^{(d)})$, i.e. $v^{(d)} = 0$. Fixing $l \ge 2$ and applying Proposition 3.2 (f), one deduces from the equations $v_l^{(d)} = 0$ for $2 \le d \le l$ that $v_{i,l-i} = 0$ for $1 \le i \le l-1$, and thus $v = 0$, ending the proof of (44).

Now, Proposition 4.1 show that $V'_d$ is stable by $\tilde{\Pi}$ and $\tilde{\Pi}'$ and, more precisely,

$$ \forall u \in L^2(\mathbb{N}_d, \mu^{(d)}), \quad \tilde{\Pi}[v(d, u)] = v(d, \Pi_d u) \quad \text{and} \quad \tilde{\Pi}'[v(d, u)] = v(d, \Pi'_d u). $$

It then follows from (44) that

$$ \|\tilde{\Pi} - \tilde{\Pi}'\| = \sup_{d \ge 2} \|\Pi_d - \Pi'_d\|_d. $$

Let us now come to the proof of (42). Let $d \ge 1$ be fixed and define

$$ g^{(d)} = \frac{d\mu^{(d+1)}}{d\mu^{(d)}} \quad \text{and} \quad h^{(d)} = \frac{d(\Pi_{d+1} - \Pi'_{d+1})}{d(\Pi_d - \Pi'_d)}, $$
i.e.

$$ \forall n \ge d+1, \quad g_n^{(d)} := \frac{(n-d)(n+d)}{2d(2d+1)} $$

and

$$ \forall i, j \ge d+1, \quad h_{i,j}^{(d)} := \begin{cases} \frac{j+d}{i+d} & \text{if } j > i, \\ \frac{j-d}{i-d} & \text{if } j < i, \\ 1 & \text{if } i = j. \end{cases} $$

For any $u \in L^2(\mathbb{N}_{d+1}, \mu^{(d+1)})$, we get

$$ \begin{aligned} \left| \frac{\langle u, (\Pi_{d+1} - \Pi'_{d+1})u \rangle_{\mu^{(d+1)}}}{\|u\|_{\mu^{(d+1)}}^2} \right| &\le \frac{\sum_{i,j \in \mathbb{N}} \mu_i^{(d)} (\Pi_d - \Pi'_d)_{i,j} g_i^{(d)} h_{i,j}^{(d)} |u_i| |u_j|}{\sum_{i \in \mathbb{N}} \mu_i^{(d)} g_i^{(d)} u_i^2} \\ &= \frac{\langle \tilde{u}, (\Pi_d - \Pi'_d) \tilde{u} \rangle_{\mu^{(d)}}}{\|u\|_{\mu^{(d)}}^2} \sup_{i,j \in \mathbb{N}_{d+1}} \sqrt{\frac{g_i^{(d)}}{g_j^{(d)}}} h_{i,j}^{(d)}, \end{aligned} $$

where

$$ \forall i \ge d, \quad \tilde{u}_i := \sqrt{g_i^{(d)}} |u_i(x)|, $$

and $\tilde{u}_i = 0$ if $i = d$. It is clear that $\tilde{u} \in L^2(\mathbb{N}_d, \mu^{(d)})$, so, taking the supremum over $u \in L^2(\mathbb{N}_{d+1}, \mu^{(d+1)})$, we have

$$ \|\Pi_{d+1} - \Pi'_{d+1}\|_{d+1} \le \sup_{i,j \in \mathbb{N}_{d+1}} \sqrt{\frac{g_i^{(d)}}{g_j^{(d)}}} h_{i,j}^{(d)} \|\Pi_d - \Pi'_d\|_d. $$

Hence, it only remains to show that

$$ \sup_{i,j \in \mathbb{N}_{d+1}} \sqrt{\frac{g_i^{(d)}}{g_j^{(d)}}} h_{i,j}^{(d)} \le 1. $$

Since $\sqrt{\frac{g_i^{(d)}}{g_i^{(d)}}} h_{i,i}^{(d)} = 1$ for all $i \ge 1$, it is sufficient to prove that

$$ \forall j > i \ge d+1, \quad \sqrt{\frac{g_i^{(d)}}{g_j^{(d)}}} h_{i,j}^{(d)} \le 1 $$
since the l.h.s. is symmetrical in $i, j$. For $j > i \ge d + 1$, we compute

$$ \sqrt{\frac{g_i^{(d)}}{g_j^{(d)}}} h_{i,j}^{(d)} = \sqrt{\frac{(i-d)(j+d)}{(j-d)(i+d)}} \le 1, $$

ending the proof of (42). A similar computation using the fact that

$$ g^{(0)} := \frac{d\mu^{(1)}}{d\mu} \quad \text{and} \quad h^{(0)} := \frac{d(\Pi_1 - \Pi'_1)}{d(\tilde{\Pi}_0 - \tilde{\Pi}'_0)} $$

are given by

$$ \forall i, j \ge 1, \quad g_i^{(0)} = 2i^2 \quad \text{and} \quad h_{i,j}^{(0)} = \frac{j}{i} $$

leads to (43).

<figure><img src="image_3.png" /></figure>

Again, let $\tilde{\Pi}^{(N)}$ be the restriction of $\tilde{\Pi}^{(N)}$ to $\mathcal{T}_N^*$, which can thus be seen as a finite $\mathcal{T}_N^* \times \mathcal{T}_N^*$ matrix. Similarly to the remark after the proof of Lemma 4.5, the spectrum of $\tilde{\Pi}^{(N)}$ is the spectrum of $\tilde{\Pi}^{(N)}$ plus the eigenvalue 0.

### 4.2.2 Spectral decomposition of infinite, compact, reversible N2dMC

The following result is an immediate consequence of Lemma 4.5 and Proposition 4.6.

**Corollary 4.7** If $\tilde{\Pi}_0$ is compact and reversible, the same is true for $\tilde{\Pi}$.

We can now extend Theorem 4.2 to the infinite compact, reversible case.

**Theorem 4.8** Assume that $Z$ has an infinite state space, i.e. $S_Z = \mathbb{Z}_+$.

(a) Theorem 4.2 (a) also holds true in this case.

(b) If $\tilde{\Pi}_0$ is compact and reversible w.r.t. the measure $(\mu_n)_{n \in \mathbb{N}}$, then, there exists an orthonormal basis of $L^2(\mathbb{N}^2, \nu)$ of right eigenvectors of $\tilde{\Pi}$ of the form (37) for $d \ge 2$, where $\nu$ is defined in (38). Moreover, any right eigenvector of $\tilde{\Pi}$ in $L^2(\mathbb{N}^2, \nu)$ is a linear combination of eigenvectors of the form (37) all corresponding to the same eigenvalue.
(c) Under the same assumptions as (b), there exists a family of right eigen-vectors of $\Pi$ of the form

$$ \{1\} \cup \left\{ P_1^{(1)}(i,j)u_{i+j}^{(1),l}, P_1^{(2)}(i,j)u_{i+j}^{(1),l} \right\}_{l \ge 1} \cup \bigcup_{d \ge 2} \left\{ P_d(i,j)u_{i+j}^{(d),l} \right\}_{l \ge 1}, \quad (45) $$

which is a basis of the vector space

$$ V := \left\{ v \in \mathbb{R}^{\mathbb{Z}_+^2} : v_{i,j} = a + \frac{i}{i+j} v_{i+j}^{(1)} + \frac{j}{i+j} v_{i+j}^{(2)} + v_{i,j}^{(3)}, \forall i, j \in \mathbb{Z}_+, \right. \\ \left. \text{with } a \in \mathbb{R}, v^{(1)}, v^{(2)} \in \mathbb{L}^2(\mathbb{N}, \mu) \text{ and } v^{(3)} \in \mathbb{L}^2(\mathbb{N}^2, \nu) \right\} \quad (46) $$

in the sense that, for all $v \in V$, there exist unique sequences $\{\alpha_l\}_{l \ge 1}$, $\{\beta_l\}_{l \ge 1}$, and $\{\gamma_{dl}\}_{d \ge 2, l \ge 1}$ and a unique $a \in \mathbb{R}$ such that

$$ v_{i,j} = a + \sum_{l \ge 1} \alpha_l i u_{i+j}^{(1),l} + \sum_{l \ge 1} \beta_l j u_{i+j}^{(1),l} + \sum_{d \ge 2} \sum_{l \ge 1} \gamma_{dl} P_d(i,j) u_{i+j}^{(d),l}, \quad (47) $$

where the series $\sum_l \alpha_l k u_k^{(1),l}$ and $\sum_l \beta_l k u_k^{(1),l}$ both converge for $\| \cdot \|_{\mu}$ and $\sum_{d,l} \gamma_{dl} P_d(i,j) u_{i+j}^{(d),l}$ converges for $\| \cdot \|_{\nu}$.

**Example 1** Assume that $Z_n = X_n + Y_n$ is a birth and death process, i.e. that the matrix $\Pi_0$ is tridiagonal. Assume moreover that all the entries just above or below the diagonal are positive (except of course for $p_{0,1}$, which is 0 since $p_{0,0} = 1$). It is well-known in this case that there always exists a reversible measure $\mu$ for $\tilde{\Pi}_0$. A well-known sufficient condition for the compactness of $\tilde{\Pi}_0$ is the case where this operator is Hilbert-Schmidt, which translates in our reversible, discrete case as

$$ \sum_{i,j \in \mathbb{N}} p_{i,j} p_{j,i} < \infty. \quad (48) $$

For a birth and death process, letting $p_k$ (resp. $q_k$) denote the birth (resp. death) probability in state $k$ and $r_k := 1 - p_k - q_k \ge 0$, this gives

$$ \sum_{i \ge 0} (r_i^2 + p_i q_{i+1}) < \infty. $$

As a side remark, note that $\tilde{\Pi}$ is not necessarily Hilbert-Schmidt when $\tilde{\Pi}_0$ is, as the condition $\sum_{(i,j),(k,l)} \pi_{(i,j),(k,l)} \pi_{(k,l),(i,j)} < \infty$ is not equivalent to (48).
**Proof** Point (a) can be proved exactly the same way as Theorem 4.2 (a).

The fact that compact selfadjoint operators admit an orthonormal basis of eigenvectors is classical. To prove (b), we only have to check that all these eigenvectors can be taken of the form (37) for $d \ge 2$. This easily follows from the fact that $\Pi_d$ is compact selfadjoint in $L^2(\mathbb{N}_d, \mu^{(d)})$ for all $d \ge 2$, from (44) and from Point (a).

The proof of (c) is similar to the proof of Theorem 2.4 (b). Fix $v \in V$ and define $a$, $v^{(1)}$, $v^{(2)}$ and $v^{(3)}$ as in (46). Since $\mu_k^{(1)} = k^2\mu_k/2$, $(v_k^{(1)}/k) \in L^2(\mathbb{N}, \mu^{(1)})$ and there exists $(\alpha_l)_{l \ge 1}$ such that

$$v_k^{(1)} = k \sum_{l \ge 1} \alpha_l u_k^{(1),l} = \sum_{l \ge 1} \alpha_l k u_k^{(1),l}, \quad \forall k \in S_Z^*, \quad (49)$$

where the convergence of the first series holds for $\| \cdot \|_{\mu^{(1)}}$, and thus of the second series for $\| \cdot \|_{\mu}$. A similar decomposition for $v^{(2)}$ and the use of (b) for $v^{(3)}$ complete the proof of (47).

It only remains to observe that, for all $v \in V$, the equation

$$v_{i,j} = a + \frac{i}{i+j} v_{i+j}^{(1)} + \frac{j}{i+j} v_{i+j}^{(2)} + v_{i,j}^{(3)}$$

uniquely characterizes $a \in \mathbb{R}$, $v^{(1)}$, $v^{(2)} \in L^2(\mathbb{N}, \mu)$ and $v^{(3)} \in L^2(\mathbb{N}^2, \nu)$. Indeed, since $v_{0,j}^{(3)} = v_{i,0}^{(3)} = 0$, one must have $a = v_{0,0}$, $v_i^{(1)} = v_{i,0} - a$ and $v_j^{(2)} = v_{0,j} - a$. $\square$

# 5 On Dirichlet eigenvalues in the complement of triangular subdomains

In this section, we consider the same model as in the previous section, and we assume either that $S = \mathcal{T}_N$ is finite or that $S = \mathbb{Z}_+^2$ and the restriction $\tilde{\Pi}_0$ of $\Pi_0$ to $\mathbb{N}$ is compact and reversible w.r.t. some measure $\mu$. We recall that $\mathcal{T}_k = \{(i,j) \in \mathbb{Z}_+^2 : i+j \le k\}$ and $\mathcal{T}_k^* = \mathcal{T}_k \cap \mathbb{N}^2$ and we define $S_k^* := S^* \setminus \mathcal{T}_{k-1}^*$ for all $k \ge 2$. Note that $S_2^* = S^*$. We also define $S_1^* := S \setminus \{(0) \times \mathbb{Z}_+\}$. Finally, for $k \ge 0$, we call $\tilde{\Pi}_k$ the restriction of the matrix $\Pi_0$ to $\{i \in S_Z : i > k\}$, and for $k \ge 1$, $\hat{\Pi}_k$ the restriction of the matrix $\Pi$ to $S_k^*$. Note that this notation is consistent with the previous definition of $\tilde{\Pi}_0$ and that $\hat{\Pi}_2 = \hat{\Pi}$. Again, all the notations of this section are gathered for reference in Appendix A.
## 5.1 The case of finite state space

Let us first assume that $\mathcal{S} = \mathcal{T}_N$ for some $N \ge 1$.

For $1 \le k \le N$, the Dirichlet eigenvalue problem for the matrix $\Pi$ in the set $S_k^*$ consists in finding $\theta \in \mathbb{C}$ and $v$ in $\mathbb{R}^{\mathcal{T}_N}$, such that

$$ \left\{ \begin{array}{ll} (\Pi v)_{(i,j)} = \theta v_{(i,j)} & \forall (i,j) \in S_k^*, \\ v_{(i,j)} = 0 & \forall (i,j) \in S \setminus S_k^*. \end{array} \right. \qquad (50) $$

This is equivalent to finding a right eigenvector of $\hat{\Pi}_k$ and extending this vector by 0 to indices in $\mathcal{S} \setminus S_k^*$. For all $k \ge 1$, we define $\theta_k^D$ as the supremum of the moduli of all Dirichlet eigenvalues in $S_k^*$. By Perron-Fröbenius' theory, $\theta_k^D$ is a Dirichlet eigenvalue in $S_k^*$.

For all $d \ge 0$, we also define $\theta^{(d)}$ as the supremum of the moduli of all eigenvalues of $\Pi$ corresponding to right eigenvectors of the form $P_d(i, j)u_{i+j}$. Again, by Theorem 4.2 (a) and Perron-Fröbenius' theory, $\theta^{(d)}$ is an eigenvector of $\Pi$ of the form $P_d(i, j)u_{i+j}$. Note that, for all $d \ge 2$, because of Proposition 3.2 (c) and (e), any right eigenvector of $\Pi$ of the form $P_d(i, j)u_{i+j}$ is a Dirichlet eigenvector of $\Pi$ in $S_{d'}$ for all $1 \le d' \le d$. In particular, $\theta_{d'}^D \ge \theta^{(d)}$. The next result gives other inequalities concerning these eigenvalues.

**Theorem 5.1** Assume that $\mathcal{S} = \mathcal{T}_N$.

(a) Then,

$$ \begin{array}{ccccccc} \theta_1^D & \ge & \theta_2^D & \ge & \theta_3^D & \ge \dots & \theta_{N-1}^D & \ge \theta_N^D \\ \parallel & & \parallel & & \vee| & & \vee| & \parallel \\ 1 = \theta^{(0)} & \ge & \theta^{(1)} & \ge & \theta^{(2)} & \ge & \theta^{(3)} & \ge \dots \ge \theta^{(N-1)} & \ge \theta^{(N)} = p_{N,N}. \end{array} $$

(b) If $\tilde{\Pi}_{k-1}$ is irreducible for some $1 \le k \le N-1$, then

$$ \theta_k^D > \theta_{k+1}^D, \quad \theta^{(k)} > \theta^{(k+1)} \qquad (51) $$

and

$$ \theta_k^D > \theta^{(k)}, \quad \text{if } k \ge 3. \qquad (52) $$

If $\tilde{\Pi}_0$ is irreducible and $p_{i,0} > 0$ for some $1 \le i \le N$, then $\theta^{(1)} < 1$.
**Proof** Since the matrix $\Pi_0$ is stochastic, it is clear from Theorem 4.2 (a) that $\theta^{(0)} = 1$. By (36), $p_{n,m}^{(1)} = \frac{m}{n}p_{n,m}$ for all $1 \le n, m \le N$ and thus $\tilde{\Pi}_0 u = \theta u$ if and only if $\Pi_1 v = \theta v$, where $v_n = u_n/n$. Hence, the largest eigenvalue of $\tilde{\Pi}_0$ is $\theta^{(1)}$. Since

$$ \Pi_0 = \begin{pmatrix} 1 & 0 \\ v & \tilde{\Pi}_0 \end{pmatrix} $$

for some (column) vector $v \in \mathbb{R}_+^N$, it is clear that $\theta^{(1)} \le \theta^{(0)}$. Ordering conveniently the states in $S_1^*$, the matrix $\hat{\Pi}_1$ has the form

$$ \hat{\Pi}_1 = \begin{pmatrix} \tilde{\Pi}_0 & 0 \\ Q & \hat{\Pi}_2 \end{pmatrix} \qquad (53) $$

for some rectangular nonnegative matrix $Q$, since the set $\{0, 1, \dots, N\} \times \{0\}$ is absorbing for the Markov chain. Again, since this matrix is block triangular, we have $\theta_1^D = \max\{\theta^{(1)}, \theta_2^D\}$. Since $\hat{\Pi}_2 = \tilde{\Pi}$, Theorem 4.2 (b) shows that

$$ \theta_2^D = \sup_{k \ge 2} \theta^{(k)}. $$

Since in addition $\hat{\Pi}_N = p_{N,N}$Id and $\Pi_N = p_{N,N}$, Theorem 5.1 (a) will hold true if we prove that the sequences $(\theta^{(k)})_{1 \le k \le N}$ and $(\theta_k^D)_{2 \le k \le N}$ are both nonincreasing.

By Perron-Fröbenius' characterization of the spectral radius of nonnegative matrices (cf. e.g. [12]), for all $1 \le k \le N$,

$$ \theta^{(k)} = \sup_{u \in \mathbb{R}^{N_k}, u \ge 0, u \ne 0} \inf_{i \in \mathbb{N}_k} \frac{\sum_{j \in \mathbb{N}_k} p_{i,j}^{(k)} u_j}{u_i}, \qquad (54) $$

where, by convention, the fraction in the r.h.s. is $+\infty$ if $u_i = 0$. Using the notation of the proof of Proposition 4.6, for all $1 \le k \le N-1$ and $u \in \mathbb{R}_+^{N_{k+1}} \setminus \{0\}$, we have

$$ \begin{aligned} \inf_{i \in \mathbb{N}_{k+1}} \frac{\sum_{j \in \mathbb{N}_{k+1}} p_{i,j}^{(k+1)} u_j}{u_i} &= \inf_{i \in \mathbb{N}_{k+1}} \frac{\sum_{j \in \mathbb{N}_{k+1}} p_{i,j}^{(k)} \tilde{u}_j \sqrt{g_i^{(k)}/g_j^{(k)}} h_{i,j}^{(k)}}{\tilde{u}_i} \\ &\le \inf_{i \in \mathbb{N}_k} \frac{\sum_{j \in \mathbb{N}_k} p_{i,j}^{(k)} \tilde{u}_j}{\tilde{u}_i}. \end{aligned} $$
Taking the supremum over $u \in \mathbb{R}_+^{\mathbb{N}_{k+1}} \setminus \{0\}$ yields

$$ \theta^{(k+1)} \le \theta^{(k)}. \qquad (55) $$

For all $k \ge 2$, the Dirichlet eigenvectors in $S_k^*$ belong to the vector space

$$ \mathcal{U}_k := \{v \in \mathbb{R}^{\mathcal{T}_N} : v = 0 \text{ on } \mathcal{T}_N \setminus S_k^*\}. $$

By Perron-Fröbenius' theory again, for all $k \ge 2$,

$$ \theta_k^D = \sup_{w \in \mathcal{U}_k \setminus \{0\}, w \ge 0} \inf_{(i,j) \in S_k^*} \frac{\sum_{(k,l) \in S_k^*} \pi_{(i,j),(k,l)} w_{(k,l)}}{w_{(i,j)}}. \qquad (56) $$

Since $S_k^* \subset S_{k-1}^*$, for all $2 \le k \le N-1$,

$$ \theta_{k+1}^D = \sup_{w \in \mathcal{U}_{k+1} \setminus \{0\}, w \ge 0} \inf_{(i,j) \in S_k^*} \frac{\sum_{(k,l) \in S_k^*} \pi_{(i,j),(k,l)} w_{(k,l)}}{w_{(i,j)}} \le \theta_k^D. \qquad (57) $$

This ends the proof of (a).

In the case where $\tilde{\Pi}_{k-1}$ is irreducible for some $2 \le k \le N-1$, it is clear that $\tilde{\Pi}_k$ and $\Pi_k$ are both irreducible. Then, by Perron-Fröbenius' theory, $\theta_k^D$ (resp. $\theta^{(k)}$) is an eigenvalue of $\tilde{\Pi}_k$ (resp. $\Pi_k$) with multiplicity one, and the corresponding nonnegative eigenvector has all its entries positive. In addition, $\theta_k^D$ (resp. $\theta^{(k)}$) is the only eigenvalue of $\tilde{\Pi}_k$ (resp. $\Pi_k$) corresponding to a positive eigenvector. In particular, the supremum in (56) (resp. (54)) is attained only at vectors $w \in V^{(k-1)}$ having all coordinates corresponding to states in $S_{k-1}^*$ positive (resp. at vectors $u \in (0, \infty)^{\mathbb{N}_k}$). Hence, the inequalities in (57) and (55) are strict.

In the case where $\tilde{\Pi}_0$ is irreducible, the same argument shows that $\theta^{(1)} > \theta^{(2)}$, and since $\theta^{(1)} = \theta_1^D$ and $\theta^{(2)} = \theta_2^D$, we also have $\theta_1^D > \theta_2^D$. This ends the proof of (51).

In the case where $\tilde{\Pi}_{k-1}$ is irreducible for $3 \le k \le N-1$, let $(u_i)_{i \in \mathbb{N}_k}$ be a positive right eigenvector of $\Pi_k$. Then the vector $P_k(i, j)u_{i+j}$ belongs to $\mathcal{U}_k$ and its restriction to $S_k^*$ is a right eigenvector of $\tilde{\Pi}_k$. However, its has positive and negative coordinates by Proposition 3.2 (g). Therefore (52) is proved.

Finally, if $\tilde{\Pi}_0$ is irreducible and $p_{i,0} > 0$ for some $1 \le i \le N$, then the absorbing state 0 is accessible by the Markov chain from any initial state (possibly after several steps). It is then standard to prove that there exists $n$ such that the sums of the entries of each line of $(\tilde{\Pi}_0)^n$ is strictly less than 1. This proves that $(\tilde{\Pi}_0)^n$ cannot have 1 as eigenvalue, and thus $\theta^{(1)} < 1$. $\square$
## 5.2 The case of infinite state space

Our goal here is to extend the previous result to the case where $S = \mathbb{Z}_+^2$ and $\tilde{\Pi}_0$ is compact reversible. So let us assume that $\tilde{\Pi}_0$ is compact and is reversible w.r.t. some measure $\mu$.

In the case of infinite state space, when $k \ge 2$ the Dirichlet eigenvalue problem for $\Pi$ in $S_k^*$ consists in finding $\theta \in \mathbb{C}$ and $v \in \mathbb{L}^2(\mathbb{Z}_+^2, \nu)$ (where the measure $\nu$ is extended by convention by zero on $\mathbb{Z}_+ \times \{0\} \cup \{0\} \times \mathbb{Z}_+$) satisfying (50). Defining the vector space where Dirichlet eigenvectors are to be found

$$ \mathcal{U}_k := \{v \in \mathbb{L}^2(\mathbb{Z}_+^2, \nu) : v = 0 \text{ on } \mathbb{Z}_+^2 \setminus S_k^*\}, $$

the supremum $\theta_k^D$ of the moduli of all Dirichlet eigenvalues in $S_k^*$ is given by

$$ \theta_k^D = \sup_{u \in \mathcal{U}_k \setminus \{0\}} \frac{|\langle u, \tilde{\Pi} u \rangle_\nu|}{\|u\|_\nu^2}. \qquad (58) $$

In view of Theorem 4.8 (c), the natural space to define the Dirichlet eigenvalue problem in $S_1^*$ is

$$ \mathcal{U}_1 := \left\{ v \in \mathbb{R}^{\mathbb{Z}_+^2} : v = 0 \text{ in } \{0\} \times \mathbb{Z}_+ \text{ and } v_{i,j} = \frac{i}{i+j} v_{i+j}^{(1)} + v_{i,j}^{(3)} \quad \forall (i,j) \in S_1^*, \right. \\ \left. \text{where } v^{(1)} \in \mathbb{L}^2(\mathbb{N}, \mu) \text{ and } v^{(3)} \in \mathbb{L}^2(\mathbb{N}^2, \nu) \right\}, $$

equipped with the norm $\| \cdot \|_{\mathcal{U}_1}$, where $\|v\|_{\mathcal{U}_1}^2 = \|v^{(1)}\|_{\mu}^2 + \|v^{(3)}\|_{\nu}^2$ (this norm is well-defined since $v^{(1)}$ and $v^{(3)}$ are uniquely defined from $v \in \mathcal{U}_1$). Then, the Dirichlet eigenvalue problem in $S_1^*$ consists in finding $\theta \in \mathbb{C}$ and $v \in \mathcal{U}_1$ satisfying (50). We also define $\theta_1^D$ as the supremum of the moduli of all Dirichlet eigenvalues in $S_1^*$.

For all $d \ge 1$, we also define $\theta^{(d)}$ as the supremum of the moduli of all eigenvalues of $\Pi$ corresponding to right eigenvectors of the form $P_d(i, j)u_{i+j}$ with $u \in \mathbb{L}^2(\mathbb{N}_d, \mu^{(d)})$. By Theorem 4.8 (a),

$$ \theta^{(d)} = \sup_{u \in \mathbb{L}^2(\mathbb{N}_d, \mu^{(d)}) \setminus \{0\}} \frac{|\langle u, \Pi_d u \rangle_{\mu^{(d)}}|}{\|u\|_{\mu^{(d)}}^2}, \quad \forall d \ge 0. $$

In addition, using the notation $\mathcal{V}_d'$ defined in the proof of Proposition 4.6, it follows from (44) that

$$ \theta^{(d)} = \sup_{u \in \mathcal{V}_d'} \frac{|\langle u, \tilde{\Pi} u \rangle_\nu|}{\|u\|_\nu^2}, \quad \forall d \ge 2. $$
Comparing this with (58), we again deduce from Proposition 3.2 (c) and (e) that $\theta_{d'}^D \ge \theta^{(d)}$ for all $2 \le d' \le d$.

Finally, since the matrix $\Pi_0$ is not reversible, we need to define $\theta^{(0)}$ in a slightly different way: $\theta^{(0)}$ is the supremum of the moduli of all eigenvalues of $\Pi$ corresponding to right eigenvectors of the form $P_d(i, j)u_{i+j}$ with $u = a\mathbf{1} + v$, where $a \in \mathbb{R}$, $\mathbf{1}$ is the vector of $\mathbb{R}^{\mathbb{Z}_+}$ with all coordinates equal to 1, $v \in \mathbb{L}^2(\mathbb{N}, \mu)$ with the convention $v_0 = 0$.

**Theorem 5.2** Assume that $S = \mathbb{Z}_+^2$ and that $\tilde{\Pi}_0$ is compact and reversible w.r.t. a positive measure $\mu$. Then, for all $d \ge 1$, $\theta_d^D$ is a Dirichlet eigenvalue of $\Pi$ in the set $S_d^*$ and $\theta^{(d)}$ is a right eigenvalue of $\Pi$ for an eigenvector of the form $P_d(i, j)u_{i+j}$ with $u \in \mathbb{L}^2(\mathbb{N}_d, \mu^{(d)})$. In addition,

$$ \begin{array}{ccccccc} \theta_1^D & \ge & \theta_2^D & \ge & \theta_3^D & \ge & \theta_4^D & \ge \dots \\ \parallel & & \parallel & & \vee| & & \vee| & \\ 1 = \theta^{(0)} & \ge & \theta^{(1)} & \ge & \theta^{(2)} & \ge & \theta^{(3)} & \ge \dots \end{array} $$

**Proof** For all $k \ge 1$, the fact that $\theta^{(k)}$ is an eigenvalue of $\Pi$ for a right eigenvector of the form $P_k(i, j)u_{i+j}$ with $u \in \mathbb{L}^2(\mathbb{N}_k, \mu^{(k)})$ follows from Theorem 4.8 (a) and from the fact that $\Pi_k$ is compact reversible. Indeed, $\|\Pi_k\|_k$ is an eigenvalue of the compact reversible matrix $\Pi_k$ (the corresponding eigenvector can be obtained for example as the limit of nonnegative eigenvectors of finite dimensional projections of $\Pi_k$). The result follows since, necessarily, $|\theta^{(k)}| \le \|\Pi_k\|_k$.

For all $k \ge 2$, we define $\Pi^{(k)}$ as the matrix whose restriction to $S_k^*$ is $\hat{\Pi}_k$ and with all other coordinates equal to zero. As the projection of a compact operator, this matrix is compact. Since it is trivially reversible for the measure $\nu$, Theorem 4.8 (b) applies to $\Pi^{(k)}$. Then, Proposition 4.6 tells us that $\theta_k^D = \|\Pi_2^{(k)}\|_2$, where $\Pi_2^{(k)}$ is defined from $\Pi^{(k)}$ as $\Pi_2$ has been defined from $\Pi$. Therefore, the fact that $\theta_k^D$ is a Dirichlet eigenvalue for $\Pi$ in $S_k^*$ can be deduced exactly as above.

Recall the definition of $\tilde{\Pi}^{(N)}$ in Section 4.2.1. For any $N \in \mathbb{N}$, replacing $\tilde{\Pi}$ by $\tilde{\Pi}^{(N)}$, we define similarly as above the quantities $\theta_k^{D,(N)}$ and $\theta^{(k,N)}$. Due to the remark after the proof of Proposition 4.6, we are again brought back to the finite framework. The following result is an immediate consequence of Lemma 4.5 and Proposition 4.6.
**Lemma 5.3** We have

$$ \forall k \ge 2, \quad \lim_{N \to \infty} \theta_k^{D,(N)} = \theta_k^D $$

$$ \forall k \ge 1, \quad \lim_{N \to \infty} \theta^{(k,N)} = \theta^{(k)} $$

From this lemma and Theorem 5.1 immediately follow all the inequalities in Theorem 5.2 that concern $\theta_k^D$ for $k \ge 2$ and $\theta^{(k)}$ for $k \ge 1$.

As in the finite case, we easily deduce from the facts that $p_{i,j}^{(1)} = \frac{i}{i}p_{i,j}$ and $\mu_i^{(1)} = 2i^2\mu_i$ that $\tilde{\Pi}_0u = \theta u$ with $u \in L^2(\mathbb{N}, \mu)$ iff $\Pi_1v = \theta v$ with $v_i = u_i/i$ and $v \in L^2(\mathbb{N}, \mu^{(1)})$. Therefore,

$$ \theta^{(1)} = \sup_{u \in L^2(\mathbb{N}, \mu) \setminus \{0\}} \frac{\langle u, \tilde{\Pi}_0 \rangle_{\mu}}{\|u\|_{\mu}}. $$

Since for all $a \in \mathbb{R}$ and $v \in L^2(\mathbb{N}, \mu)$, $\Pi_0(a\mathbf{1} + v) = a\mathbf{1} + \tilde{\Pi}_0v$, we deduce that $\theta^{(0)} = \sup\{1; \theta^{(1)}\}$. Since $\tilde{\Pi}_0$ is substochastic and reversible w.r.t. $\mu$, we have for all $u \in L^2(\mathbb{N}, \mu)$

$$ \frac{\langle u, \tilde{\Pi}_0 \rangle_{\mu}}{\|u\|_{\mu}} = \frac{\sum_{i,j \ge 1} \mu_i p_{i,j} u_i u_j}{\sum_{i \ge 1} u_i^2 \mu_i} \le \frac{\sqrt{\sum_{i,j} u_i^2 \mu_i p_{i,j}} \sqrt{\sum_{i,j} u_j^2 \mu_i p_{i,j}}}{\sum_{i \ge 1} u_i^2 \mu_i} \le 1. $$

This yields

$$ \theta^{(0)} = 1 \ge \theta^{(1)}. $$

In order to complete the proof, it only remains to check that $\theta_1^D = \theta^{(1)}$ and that $\theta_1^D$ is a Dirichlet eigenvalue in $S_1^*$. As in the finite case, $\tilde{\Pi}_1$ has the block-triangular form (53). Therefore, we obviously have $\theta_1^D \ge \theta_2^D$. In addition, any Dirichlet eigenvalue in $S_1^*$ which corresponds to an eigenvector in $U_1$ which is nonzero on the set of indices $\mathbb{N} \times \{0\}$, must be an eigenvalue of $\tilde{\Pi}_0$ corresponding to a right eigenvector in $L^2(\mathbb{N}, \mu)$. Now, if $u \in L^2(\mathbb{N}, \mu)$ satisfies $\tilde{\Pi}_0u = \theta^{(1)}u$, then $\Pi_1v = \theta^{(1)}v$ with $v_i = u_i/i$ and it follows from Theorem 4.8 (a) that the vector $\frac{i}{i+j}u_{i+j}$ is a right eigenvector of $\tilde{\Pi}_1$. Since this vector obviously belongs to $U_1$, we obtain that $\theta_1^D = \theta^{(1)}$ and that $\theta_1^D$ is a Dirichlet eigenvalue in $S_1^*$. $\square$
# 6 Application to quasi-stationarity in nearly neutral finite absorbed two-dimensional Markov chains

In this section, we restrict ourselves to the finite state space case for simplicity: let $S = T_N$ for some $N \in \mathbb{N}$. The first coordinate will be referred to as **type 1** and the second coordinate as **type 2**. Recall that the sets $\mathbb{Z}_+ \times \{0\}$, $\{0\} \times \mathbb{Z}_+$ and $\{(0,0)\}$ are absorbing for the N2dMC considered above, which means that each sub-population in the model can go extinct. This means that the transition matrix $\Pi$ has the form

$$ \Pi = \begin{pmatrix} 1 & 0 \\ r & Q \end{pmatrix}, \qquad (59) $$

after ordering the states as (0,0) first. Ordering the states in $S \setminus \{0\}$ as $\{1, \dots, N\} \times \{0\}$ first and $\{0\} \times \{1, \dots, N\}$ second, the matrix $Q$ has the form

$$ Q = \begin{pmatrix} Q_1 & 0 & 0 \\ 0 & Q_2 & 0 \\ R_1 & R_2 & Q_3 \end{pmatrix} \qquad (60) $$

where $Q_i$ ($1 \le i \le 3$) are square matrices and $R_i$ ($1 \le i \le 2$) rectangular matrices.

In this section, we study the problem of quasi-stationary distributions (QSD) and quasi-limiting distributions (QLD, see the introduction) for Markov processes, not necessarily neutral, whose transition matrix has the form (59–60). The classical case [7] for such a study is the case when $Q$ is irreducible, which does not hold here. A general result is proved in Subsection 6.1. Our results of Section 5 are then applied in Subsection 6.2 to study the quasi-limiting distribution of nearly neutral two-dimensional Markov chains.

## 6.1 Yaglom limit for general absorbing two-dimensional Markov chains

Let $(X_n, Y_n, n \ge 0)$ be a Markov chain on $S = T_N$, with transition matrix of the form (59–60). We do not assume that this process is neutral. We call such processes A2dMC for “absorbed two-dimensional Markov chains”.

Under the assumption that the matrices $Q_1, Q_2$ and $Q_3$ are irreducible and aperiodic, Perron-Frobenius' theory ensures the existence of a unique
eigenvalue $\theta_i$ with maximal modulus for $Q_i$, which is real and has multiplicity one. Moreover, $Q_i$ admits unique (strictly) positive right and left eigenvectors for $\theta_i$, $u_i$ and $v_i$ respectively, normalised as $v_i \mathbf{1} = 1$ and $v_i u_i = 1$, where we use the convention that $u_i$ is a row vector and $v_i$ is a column vector, and where $\mathbf{1}$ denotes the column vector with adequate number of entries, all equal to 1. In the following result, we use the classical identification of column vectors and measures: for example, $v_1 = ((v_1)_i)_{(i,0)\in S\setminus\{(0,0)\}}$ is identified with $\sum_{(i,0)\in S\setminus\{(0,0)\}} (v_1)_i \delta_i$, and $v_3 = ((v_3)_{(i,j)})_{(i,j)\in S^*}$ is identified with $\sum_{(i,j)\in S\cap\mathbb{N}^2} (v_3)_{(i,j)} \delta_{(i,j)}$.

With this notation, using the results of Darroch and Seneta [7], $v_1 \otimes \delta_0$ and $\delta_0 \otimes v_2$ are trivial QSDs for the Markov chain $(X, Y)$.

**Theorem 6.1** Assume that the matrices $Q_1, Q_2$ and $Q_3$ are irreducible and aperiodic, $R_1 \neq 0$ and $R_2 \neq 0$. Then, for any $i \ge 1$ such that $(i, 0) \in S$,

$$ \lim_{n \to +\infty} \mathcal{L}_{(i,0)}[(X_n, Y_n) | (X_n, Y_n) \neq (0,0)] = v_1 \otimes \delta_0, \quad (61) $$

and similarly for the initial state $(0, i) \in S$, where $\mathcal{L}_{(i,j)}$ denotes the law of the Markov chain $(X_n, Y_n)_{n \ge 0}$ with initial condition $(i, j)$.

Moreover, for any $(i, j) \in S^*$,

$$ \lim_{n \to +\infty} \mathcal{L}_{(i,j)}[(X_n, Y_n) | (X_n, Y_n) \neq (0,0)] = \begin{cases} v_1 \otimes \delta_0 & \text{if } \theta_1 \ge \theta_3 \text{ and } \theta_1 > \theta_2, \\ \delta_0 \otimes v_2 & \text{if } \theta_2 \ge \theta_3 \text{ and } \theta_2 > \theta_1, \\ \frac{v_3 + w_1 \otimes \delta_0 + \delta_0 \otimes w_2}{1 + w_1 \mathbf{1} + w_2 \mathbf{1}} & \text{if } \theta_3 > \theta_1, \theta_2, \\ p_{i,j} v_1 \otimes \delta_0 + (1 - p_{i,j}) \delta_0 \otimes v_2 & \text{if } \theta_1 = \theta_2 > \theta_3, \\ q v_1 \otimes \delta_0 + (1 - q) \delta_0 \otimes v_2 & \text{if } \theta_1 = \theta_2 = \theta_3. \end{cases} \quad (62) $$

where

$$ w_i = v_3 R_i (\theta_3 I - Q_i)^{-1}, \quad i = 1, 2, \quad (63) $$

$$ p_{i,j} = \frac{\delta_{(i,j)} (\theta_1 I - Q_3)^{-1} R_1 u_1}{\delta_{(i,j)} (\theta_1 I - Q_3)^{-1} (R_1 u_1 + R_2 u_2)} \quad (64) $$

$$ \text{and} \quad q = \frac{v_3 R_1 u_1}{v_3 (R_1 u_1 + R_2 u_2)}. \quad (65) $$
To give an interpretation of this result, we will say that there is **extinction** of type *i* conditionally on non-extinction if the QLD (62) gives mass 0 to all states with positive *i*-th coordinate. Conversely, we say that there is **coexistence** conditionally on non-extinction if the QLD (62) gives positive mass to the set $S^*$. We also say that type 1 is **stronger** than type 2 (or type 2 is **weaker** than type 1) if $\theta_1 > \theta_2$, and conversely if $\theta_2 > \theta_1$.

Theorem 6.1 says that the limit behaviour of the population conditionally on non-extinction essentially depends on whether the largest eigenvalue of $Q$ is $\theta_1$, $\theta_2$ or $\theta_3$. If either $\theta_1 > \theta_2$ and $\theta_1 \ge \theta_3$ or $\theta_2 > \theta_1$ and $\theta_2 \ge \theta_3$, the QLD is the same as if there were no individual of the weaker type in the initial population, and there is extinction of the weaker type conditionally on non-extinction. If $\theta_3 > \theta_1$, $\theta_2$, there is **coexistence** of both types conditionally on non-extinction. Finally, when $\theta_1 = \theta_2$, both types can survive under the QLD, so none of the types go extinct (according to the previous terminology), but there is no coexistence, as one (random) type eventually goes extinct. Observe also that the case $\theta_1 = \theta_2 > \theta_3$ is the only one where the QLD depends on the initial condition.

Note also that, in the case where $\theta_3 \le \max\{\theta_1, \theta_2\}$ and $\theta_1 \ne \theta_2$, the QLD does not depend on any further information about the matrix $Q_3$. In other words, if one knows a priori that $\theta_3 \le \max\{\theta_1, \theta_2\}$ and $\theta_1 \ne \theta_2$, the precise transition probabilities of the Markov chain from any state in $S^*$ have no influence on the QLD. The QLD is only determined by the monotype chains of types 1 and 2.

Our next result says that, for any values of $\theta_1$, $\theta_2$ and $\theta_3$, all the QSDs of the Markov chain are those given in the r.h.s. of (62), when they exist and are nonnegative.

**Proposition 6.2** Under the same assumptions and notation as in Theorem 6.1, the set of QSDs of the Markov chain is composed of the probability measures $p v_1 \otimes \delta_0 + (1-p) \delta_0 \otimes v_2$ for all $p \in [0, 1]$, with the additional QSD

$$ \frac{v_3 + w_1 \otimes \delta_0 + \delta_0 \otimes w_2}{1 + w_1 \mathbf{1} + w_2 \mathbf{1}} \qquad (66) $$

in the case where $\theta_3 > \max\{\theta_1, \theta_2\}$.

**Proof** The fact that all the QSDs giving no mass to the set $S^*$ are of the form $p v_1 \otimes \delta_0 + (1-p) \delta_0 \otimes v_2$ for some $p \in [0, 1]$ is an immediate consequence of the facts that the sets $\{1, \dots, N\} \times \{0\}$ and $\{0\} \times \{1, \dots, N\}$ do not
communicate and the only QSD of an irreducible and aperiodic Markov chain on a finite set is given by the only positive normalized left eigenvector of the transition matrix of the Markov chain (cf. [7]).

Assume now that $\mu$ is a QSD for the Markov chain $(X_n, Y_n)_{n \ge 0}$ such that $\mu(S^*) > 0$, and write $\mu = (\mu_1, \mu_2, \mu_3)$, where $\mu_1$ (resp. $\mu_2$, resp. $\mu_3$) is the restriction of $\mu$ to the set $\{1, \dots, N\} \times \{0\}$ (resp. $\{0\} \times \{1, \dots, N\}$, resp. $S^*$). The equation $\mu Q = \theta Q$ for some $\theta > 0$, which characterizes QSDs, implies that $\mu_3$ is a nonnegative left eigenvector for $Q_3$. Thus, by the Perron-Frobenius theorem, $\mu_3 = av_3$ for some $a > 0$ and $\theta = \theta_3$. Using again the formula $\mu Q = \theta_3 Q$, one necessarily has

$$ \mu_i(\theta_3 \mathrm{Id} - Q_i) = av_3 R_i, \quad i = 1, 2. \qquad (67) $$

In the case where $\theta_3 > \max\{\theta_1, \theta_2\}$, the matrices $\theta_3 \mathrm{Id} - Q_1$ and $\theta_3 \mathrm{Id} - Q_1$ are invertible, as shown in Lemma 6.3 below. Thus $\mu$ is given by (66).

In the case where $\theta_3 \le \theta_i$ for $i = 1$ or $2$, we deduce from (67) that $(\theta_3 - \theta_i)\mu_i u_i = av_3 R_i u_i$. This is impossible since the l.h.s. of this formula is non-positive and the r.h.s. is positive as $R_i \ne 0$, $v_3 > 0$ and $u_i > 0$. $\square$

**Proof of Theorem 6.1** For all $(k, l) \in S \setminus \{(0, 0)\}$, we want to compute the limit of

$$ \mathbb{P}_{(i,j)}[(X_n, Y_n) = (k, l) | (X_n, Y_n) \ne (0, 0)] = \frac{Q_{(i,j),(k,l)}^{(n)}}{\sum_{(k', l') \ne (0,0)} Q_{(i,j),(k',l')}^{(n)}} \quad (68) $$

as $n \to +\infty$, where $Q_{(i,j),(k,l)}^{(n)}$ denotes the element of $Q^n$ on the line corresponding to state $(i, j) \in S$ and the column corresponding to state $(k, l) \in S$.

Therefore, we need to analyse the behaviour of $Q^n$ as $n \to +\infty$. We have by induction

$$ Q^n = \begin{pmatrix} Q_1^n & 0 & 0 \\ 0 & Q_2^n & 0 \\ R_1^{(n)} & R_2^{(n)} & Q_3^n \end{pmatrix} \qquad (69) $$

where

$$ R_i^{(n)} = \sum_{k=0}^{n-1} Q_3^k R_i Q_i^{n-1-k}, \quad i = 1, 2. \qquad (70) $$
By the Perron-Frobenius Theorem (see e.g. [12]),

$$Q_i^n = \theta_i^n u_i v_i + O((\theta_i \alpha)^n \mathbf{1}) \quad (71)$$

for some $\alpha < 1$ as $n \to +\infty$, where $\mathbf{1}$ denotes the square matrix of appropriate dimension, whose entries are all equal to 1. We need the following result. Its proof is postponed at the end of the subsection.

**Lemma 6.3** If $\theta_1 > \theta_3$, the matrix $\theta_1 \mathrm{Id} - Q_3$ is invertible, its inverse has positive entries and

$$R_1^{(n)} \sim \theta_1^n (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 v_1 \quad (72)$$

as $n \to +\infty$. If $\theta_1 < \theta_3$, the matrix $\theta_3 \mathrm{Id} - Q_1$ is invertible, its inverse has positive entries and

$$R_1^{(n)} \sim \theta_3^n u_3 v_3 R_1 (\theta_3 \mathrm{Id} - Q_1)^{-1} \quad (73)$$

as $n \to +\infty$. If $\theta_1 = \theta_3$, as $n \to +\infty$,

$$R_1^{(n)} \sim n \theta_1^{n-1} u_3 v_3 R_1 u_1 v_1. \quad (74)$$

Theorem 6.1 can be proved from this result as follows. Let $D$ denote the denominator of (68). If $\theta_1, \theta_2 > \theta_3$, (72) and (71) yield for all $(i, j) \in S^*$

$$D \sim \theta_1^n \delta_{(i,j)} (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 v_1 \mathbf{1} + \theta_2^n \delta_{(i,j)} (\theta_2 \mathrm{Id} - Q_3)^{-1} R_2 u_2 v_2 \mathbf{1} + \theta_3^n \delta_{(i,j)} u_3 v_3 \mathbf{1}$$

as $n \to +\infty$. In the case when $\theta_1 > \theta_2$, since $(\theta_1 \mathrm{Id} - Q_3)^{-1}$ has positive entries, we have $D \sim \theta_1^n \delta_{(i,j)} (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1$. The limit of (68) when $n \to +\infty$ then follows from (72). The case $\theta_2 > \theta_1$ is treated similarly. In the case when $\theta_1 = \theta_2$,

$$D \sim \theta_1^n \delta_{(i,j)} [(\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 + (\theta_2 \mathrm{Id} - Q_3)^{-1} R_2 u_2],$$

and the fourth line of (62) follows from Lemma 6.3.

In the case when $\theta_3 > \theta_1, \theta_2$, we obtain

$$D \sim \theta_3^n \delta_{(i,j)} u_3 v_3 [R_1 (\theta_3 \mathrm{Id} - Q_1)^{-1} \mathbf{1} + R_2 (\theta_3 \mathrm{Id} - Q_2)^{-1} \mathbf{1} + \mathbf{1}],$$

which implies the third line of (62).

Similarly, it follows from (71) and Lemma 6.3 that

$$D \sim \theta_1^n \delta_{(i,j)} (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 \quad \text{if } \theta_1 > \theta_3 > \theta_2,$$

$$D \sim n \theta_1^{n-1} \delta_{(i,j)} u_3 v_3 R_1 u_1 \quad \text{if } \theta_1 = \theta_3 > \theta_2,$$

$$D \sim n \theta_1^{n-1} \delta_{(i,j)} u_3 v_3 (R_1 u_1 + R_2 u_2) \quad \text{if } \theta_1 = \theta_2 = \theta_3.$$

The proof is easily completed in each of these cases. $\square$
**Proof of Lemma 6.3** Assume first that $\theta_1 > \theta_3$. One easily checks that

$$ (\theta_1 \mathrm{Id} - Q_3) \sum_{k=0}^{n-1} \theta_1^{-k} Q_3^k = \theta_1 \mathrm{Id} - \theta_1^{-n+1} Q_3^n. \quad (75) $$

Because of (71), the series in the previous equation converges when $n \to +\infty$. Therefore, $\theta_1 \mathrm{Id} - Q_3$ is invertible and

$$ (\theta_1 \mathrm{Id} - Q_3)^{-1} = \sum_{n \ge 0} \theta_1^{-n-1} Q_3^n, $$

which has positive entries since $Q_3$ is irreducible. Therefore, it follows from (71) and (75) that

$$ \begin{aligned} R_1^{(n)} &= \sum_{k=0}^{n-1} Q_3^k R_1 [\theta_1^{n-1-k} u_1 v_1 + O((\theta_1 \alpha)^{n-1-k} \mathbf{1})] \\ &= \theta_1^{n-1} (\theta_1 \mathrm{Id} - \theta_1^{-n+1} Q_3^n) (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 v_1 + O \left( u_3 v_3 R_1 \mathbf{1} \sum_{k=0}^{n-1} \theta_3^k (\theta_1 \alpha)^{n-1-k} \right) \\ &= \theta_1^n (\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 v_1 + O(\theta_3^n \mathbf{1}) + O((\theta_3^n + (\theta_1 \alpha)^n) \mathbf{1}), \end{aligned} $$

where we used the fact that $\alpha$ may be increased without loss of generality so that $\theta_3 \neq \theta_1 \alpha$, in which case

$$ \sum_{k=0}^{n-1} \theta_3^k (\theta_1 \alpha)^{n-1-k} = \frac{\theta_3^n - (\theta_1 \alpha)^n}{\theta_3 - \theta_1 \alpha}. $$

Since $R_1 \neq 0$ has nonnegative entries and $(\theta_1 \mathrm{Id} - Q_3)^{-1}$ and $u_1 v_1$ have positive entries, the matrix $(\theta_1 \mathrm{Id} - Q_3)^{-1} R_1 u_1 v_1$ also has positive entries, and (72) follows. The case $\theta_3 > \theta_1$ can be handled similarly.

Assume finally that $\theta_1 = \theta_3$. By (71),

$$ \begin{aligned} R_1^{(n)} &= \sum_{k=0}^{n-1} (\theta_3^k u_3 v_3 + O((\theta_3 \alpha)^k \mathbf{1})) R_1 (\theta_1^{n-k-1} u_1 v_1 + O((\theta_1 \alpha)^{n-k-1} \mathbf{1})) \\ &= n \theta_1^{n-1} u_3 v_3 R_1 u_1 v_1 + O \left( \frac{\theta_1^{n-1}}{1-\alpha} (\mathbf{1} R_1 u_1 v_1 + u_3 v_3 R_1 \mathbf{1}) \right) + O(n(\alpha \theta_1)^{n-1} \mathbf{1}), \end{aligned} $$

which ends the proof of Lemma 6.3. $\square$
## 6.2 The nearly neutral case

Since $\Pi$ in (59) is a block triangular matrix, we have

$$ \mathrm{Sp}'(\Pi) = \{1\} \cup \mathrm{Sp}'(Q_1) \cup \mathrm{Sp}'(Q_2) \cup \mathrm{Sp}'(Q_3), $$

where $\mathrm{Sp}'(A)$ denotes the spectrum of the matrix $A$, where eigenvalues are counted with their multiplicity.

In the case of a N2dMC satisfying the assumptions of Theorem 6.1, with the notation of Section 4, we have $Q_1 = Q_2 = \tilde{\Pi}_0$. By Theorem 4.2 (b) and Remark 4.3, $\{1\} \cup \mathrm{Sp}'(Q_1) \cup \mathrm{Sp}'(Q_2)$ is the set of eigenvalues corresponding to right eigenvectors of $\Pi$ of the form $P_0(i, j)u_{i+j}$ and $P_1(i, j)u_{i+j}$, counted with their multiplicity. More precisely, $\mathrm{Sp}'(Q_1)$ corresponds to eigenvectors of the form $P_1^{(1)}(i, j)u_{i+j}$, and $\mathrm{Sp}'(Q_2)$ to eigenvectors of the form $P_1^{(2)}(i, j)u_{i+j}$. In particular, $\theta_1 = \theta_2 = \theta^{(1)} = \theta_1^D$, with the notation of Theorem 5.1. Moreover, since $Q_3 = \Pi^{(2)}$, Theorem 5.1 shows that $\theta_3 = \theta^{(2)} = \theta_2^D < \theta_1 = \theta_2$ and $\mathrm{Sp}'(Q_3)$ is the set of eigenvalues corresponding to right eigenvectors of $\Pi$ of the form $P_d(i, j)u_{i+j}$ for $d \ge 2$, counted with their multiplicity.

In other words, with the terminology defined after Theorem 6.1, coexistence is impossible in the neutral case. Since the eigenvalues of $Q_1, Q_2$ and $Q_3$ depend continuously on the entries of these matrices, we deduce that coexistence is impossible in the neighborhood of neutrality:

**Corollary 6.4** Let $\Pi$ be the transition matrix of some fixed N2dMC in $\mathcal{T}_N$ such that $\tilde{\Pi}_0$ and $\tilde{\Pi}_1$ are both irreducible and there exists $i \in \{1, \dots, N\}$ such that $p_{i,0} > 0$. For any A2dMC $(X, Y)$ in $\mathcal{T}_N$ with transition matrix $\Pi'$ sufficiently close to $\Pi$, coexistence is impossible in the QLD of $(X, Y)$. Let $\theta'_1, \theta'_2, \theta'_3$ denote the eigenvalues $\theta_1, \theta_2, \theta_3$ of Theorem 6.1 corresponding to the matrix $\Pi'$. If $\theta'_1 \neq \theta'_2$, the QLD of $(X, Y)$ is the trivial QSD corresponding to the stronger type: if $\theta'_1 > \theta'_2$, the QLD of $(X, Y)$ is $v'_1 \otimes \delta_0$, where $v'_1$ is the QLD of $(X, 0)$, and if $\theta'_2 > \theta'_1$, the QLD of $(X, Y)$ is $\delta_0 \otimes v'_2$, where $v'_2$ is the QLD of $(0, Y)$.

## A Notations

We gather here all the notations used at several places in the paper. Most of these notations are introduced in Sections 4 and 5.
## A.1 General definitions

* For any measurable subset $\Gamma$ of $\mathbb{R}^d$ and any $\sigma$-finite positive measure $\mu$ on $\Gamma$, $L^2(\Gamma, \mu)$ is the set of Borel functions $f$ on $\Gamma$ defined up to $\mu$-negligible set such that $\int_{\Gamma} f^2 d\mu < +\infty$. We denote by $\langle \cdot, \cdot \rangle_{\mu}$ the canonical inner product on $L^2(\Gamma, \mu)$ and $\| \cdot \|_{\mu}$ the associated norm. In the case when $\Gamma$ is discrete, we make the usual abuse of notation to identify the measure $\mu$ and the corresponding function on $\Gamma$.

* For all $I \times I$ square matrix $M$, where $I$ is a finite or denumerable set of indices, and for all $J \subset I$, we call “restriction of the matrix $M$ to $J$” the matrix obtained from $M$ by removing all rows and columns corresponding to indices in $I \setminus J$.

## A.2 Polynomials

$H_2(X), H_3(X), \dots$ are defined in Proposition 2.2.

$$P_0(X, Y) = 1.$$

$$P_1^{(1)}(X, Y) = P_1(X, Y) = X.$$

$$P_1^{(2)}(X, Y) = Y.$$

$P_2(X, Y), P_3(X, Y), \dots$ are defined in Theorem 3.1.

## A.3 Sets

$$\mathbb{Z}_+ = \{0, 1, \dots\}.$$

$$\mathbb{N} = \{1, 2, \dots\}.$$

$$\mathbb{N}_d = \{d, d+1, \dots\}.$$

$$\mathcal{T}_N = \{(i, j) \in \mathbb{Z}_+^2 : i + j \le N\}, \text{ where } N \ge 0 \text{ is fixed below.}$$

$$\mathcal{T}_N^* = \mathcal{T}_N \cap \mathbb{N}^2.$$

$$
\begin{array}{ll}
\text{Finite case} & \text{Infinite case} \\
S_Z = \{0, 1, \dots, N\}. & S_Z = \mathbb{Z}_+ \\
S_Z^* = \{1, 2, \dots, N\}. & S_Z^* = \mathbb{N} \\
S = \mathcal{T}_N. & S = \mathbb{Z}_+^2 \\
S^* = \mathcal{T}_N^*. & S^* = \mathbb{N}^2 \\
S_k^* = \{(i, j) \in \mathbb{Z}_+^2 : k \le i + j \le N\}, & S_k^* = \{(i, j) \in \mathbb{Z}_+^2 : k \le i + j \le N\}, \text{ for all } k \ge 2 \\
\mathcal{T}_N^* = \mathcal{T}_N \cap (\mathbb{N} \times \mathbb{Z}_+) & k \ge 2 \\
S_1^* = \mathcal{T}_N \cap (\mathbb{N} \times \mathbb{Z}_+). & S_1^* = \mathbb{N} \times \mathbb{Z}_+.
\end{array}
$$
## A.4 Matrices

$\Pi_0 = (p_{n,m})_{n,m \in S_Z}$ is a stochastic matrix such that $p_{0,0} = 1$.

$\tilde{\Pi}_k$ is the restriction of $\Pi_0$ to the set of indices $S_Z \cap \mathbb{N}_{k+1}$, for all $k \ge 0$.

$\text{Pr}_N$ is the projection operator on $\mathbb{R}^N$ defined by $\text{Pr}_N(u_1, u_2, \dots) = (u_1, \dots, u_N, 0, \dots)$.

$\tilde{\Pi}_0^{(N)} = \text{Pr}_N \tilde{\Pi}_0 \text{Pr}_N$, in the infinite case (i.e. when $S_Z = \mathbb{Z}_+$).

$\Pi_0^{(N)}$ is the Markovian kernel on $\mathbb{Z}_+$ whose restriction to $\mathbb{N}$ is $\tilde{\Pi}_0^{(N)}$.

$\tilde{\Pi}_0^{(N)}$ is the restriction of $\tilde{\Pi}_0^{(N)}$ to $\{1, \dots, N\}$.

$\Pi = (\pi_{(i,j),(k,l)})(i,j),(k,l) \in S$, where

$$ \pi_{(i,j),(i+k,j+l)} = \begin{cases} \frac{\binom{i+k-1}{k} \binom{j+l-1}{l}}{\binom{i+j+k+l-1}{k+l}} p_{i+j, i+j+k+l} & \text{if } (k,l) \in \mathbb{Z}_+^2 \setminus \{0\}, \\ \frac{\binom{i}{k} \binom{j}{l}}{\binom{i+j}{k+l}} p_{i+j, i+j-k-l} & \text{if } (-k,-l) \in \mathbb{Z}_+^2, \\ 0 & \text{otherwise,} \end{cases} $$

with the convention that $\binom{i}{j} = 0$ if $i < 0$, $j < 0$ or $j > i$.

$\tilde{\Pi}$ is the restriction of $\Pi$ to $S^*$.

$\tilde{\Pi}_k$ is the restriction of the matrix $\Pi$ to $S_k^*$, for all $k \ge 1$.

$\tilde{\Pi}^{(N)}$ is constructed from $\Pi_0^{(N)}$ exactly as $\tilde{\Pi}$ is defined from $\Pi_0$.

$\tilde{\Pi}^{(N)}$ is the restriction of $\tilde{\Pi}^{(N)}$ to $T_N^*$.

$\Pi_d = (p_{n,m}^{(d)})(n,m) \in S, n,m \ge d$ for all $d \ge 0$, where for all $(n,m) \in S, n,m \ge d$,

$$ p_{n,m}^{(d)} = \begin{cases} \frac{\binom{m+d-1}{m-n}}{\binom{m-1}{m-n}} p_{n,m} & \text{if } m > n, \\ \frac{\binom{n-d}{n-m}}{\binom{n}{n-m}} p_{n,m} & \text{if } m < n, \\ p_{n,n} & \text{if } m = n. \end{cases} $$

For A2dMC (see Section 6.1)

$$ \Pi = \begin{pmatrix} 1 & 0 \\ r & Q \end{pmatrix}, \text{ where } Q = \begin{pmatrix} Q_1 & 0 & 0 \\ 0 & Q_2 & 0 \\ R_1 & R_2 & Q_3 \end{pmatrix}. $$
## A.5 Measures and vectors

$\mu = (\mu_i)_{i \in S_Z^*}$ is a reversible measure for the matrix $\tilde{\Pi}_0$.

$$ \nu = (\nu_{(i,j)})(i,j) \in S^*, \text{ where } \nu_{(i,j)} = \frac{(i+j)\mu_{i+j}}{ij}. $$

$\mu^{(d)} = (\mu_n^{(d)})_{n \in S_Z \cap \mathbb{N}_d}$ for all $d \ge 1$, where

$$ \mu_n^{(d)} = 2n \binom{n+d-1}{2d-1} \mu_n \quad \text{for all } n \in S_Z \cap \mathbb{N}_d. $$

$g^{(d)} = (g_i^{(d)})_{i \in S_Z \cap \mathbb{N}_{d+1}}$ for all $d \ge 0$, where

$$ g_i^{(0)} = 2i^2 \quad \text{for all } i \in S_Z^* \text{ and} $$

$$ g_i^{(d)} = \frac{(i-d)(i+d)}{2d(2d+1)} \quad \text{for all } d \ge 1 \text{ and } i \in S_Z \cap \mathbb{N}_{d+1}. $$

$$ h^{(d)} = (h_{i,j}^{(d)})_{(i,j) \in S, i,j \ge d+1} \quad \text{for all } d \ge 0, \text{ where} $$

$$ h_{i,j}^{(0)} = \frac{j}{i} \quad \text{for all } (i,j) \in S^* \text{ and} $$

$$ h_{i,j}^{(d)} = \begin{cases} \frac{j+d}{i+d} & \text{if } j > i, \\ \frac{j-d}{i-d} & \text{if } j < i, \\ 1 & \text{if } i = j. \end{cases} \quad \text{for all } d \ge 1 \text{ and } (i,j) \in S, i,j \ge d+1. $$

## A.6 Operator norms (in the infinite, reversible case)

$\|\cdot\|_0$ is the natural operator norm on the set of bounded operators on $L^2(\mathbb{N}, \mu)$.

$\|\cdot\|_1$ is the natural operator norm on the set of bounded operators on $L^2(\mathbb{N}^2, \nu)$.

$\|\cdot\|_d$ is the natural operator norm on the set of bounded operators on $L^2(\mathbb{N}_d, \mu^{(d)})$, for all $d \ge 1$.

## A.7 Eigenvalues

We refer to Sections 5.1 and 5.2 for precise definitions in the finite and infinite cases.

$\theta_k^D$ is the biggest Dirichlet eigenvalue of $\Pi$ in $S_k^*$, for all $k \ge 1$.

$\theta^{(d)}$ is the biggest eigenvalue of $\Pi$ corresponding to right eigenvectors of the form $P_d(i, j)u_{i+j}$, for all $d \ge 0$.

For A2dMC (see Section 6.1)
$\theta_1, \theta_2$ and $\theta_3$ are the Perron-Fröbenius eigenvalues of $Q_1, Q_2$ and $Q_3$, respectively.

## A.8 Vector spaces

$$ \mathcal{V}_d := \{v \in \mathbb{R}^S : v_{i,j} = P_d(i,j)u_{i+j} \text{ with } u \in \mathbb{R}^{S_z}\}. $$

**Finite case**

$$ \mathcal{U}_k = \{v \in \mathbb{R}^{\mathcal{T}_N} : v = 0 \text{ on } \mathcal{T}_N \setminus S_k^*\} \text{ for all } k \ge 1. $$

**Infinite, reversible case**

$$ \mathcal{U}_k = \{v \in \mathbb{L}^2(\mathbb{Z}_+^2, \nu) : v = 0 \text{ on } \mathbb{Z}_+^2 \setminus S_k^*\} \text{ for all } k \ge 2. $$

$$ \mathcal{U}_1 = \{ v \in \mathbb{R}^{\mathbb{Z}_+^2} : v_{i,j} = \frac{i}{i+j} v_{i+j}^{(1)} + v_{i,j}^{(3)}, v^{(1)} \in \mathbb{L}^2(\mathbb{N}, \mu), v^{(3)} \in \mathbb{L}^2(\mathbb{N}^2, \nu) \}. $$

$$ \mathcal{V}'_d = \{v \in \mathbb{L}^2(\mathbb{N}^2, \nu) : v_{i,j} = P_d(i,j)u_{i+j} \text{ with } u \in \mathbb{L}^2(\mathbb{N}_d, \mu^{(d)})\} \text{ for all } d \ge 2. $$

**Acknowledgments:** We are grateful to Pierre-Emmanuel Jabin for helpful discussions on the continuous case of Section 2 and to Bob Griffiths for his comments on a preliminary version of this work and for pointing out the reference [15].

Nicolas Champagnat's research was partly supported by project MANEGE 'Modèles Aléatoires en Écologie, Génétique et Évolution' 09-BLAN-0215 of ANR (French national research agency).

## References

[1] Sylvain Billiard and Viet Chi Tran. A general stochastic model for sporophytic self-incompatibility. *J. Math. Biol.*, 64(1–2):163–210, 2012.

[2] Haïm Brezis. *Analyse fonctionnelle*. Collection Mathématiques Appliquées pour la Maîtrise. [Collection of Applied Mathematics for the Master's Degree]. Masson, Paris, 1983. Théorie et applications. [Theory and applications].

[3] Patrick Cattiaux, Pierre Collet, Amaury Lambert, Servet Martínez, Sylvie Méléard, and Jaime San Martín. Quasi-stationary distributions and diffusion models in population dynamics. *Ann. Probab.*, 37(5):1926–1969, 2009.
[4] Patrick Cattiaux and Sylvie Méléard. Competitive or weak cooperative stochastic lotka-volterra systems conditioned on non-extinction. *J. Math. Biol.*, 2010. To appear.

[5] James A. Cavender. Quasi-stationary distributions of birth-and-death processes. *Adv. Appl. Probab.*, 10(3):570–586, 1978.

[6] Nicolas Champagnat and Amaury Lambert. Evolution of discrete populations and the canonical diffusion of adaptive dynamics. *Ann. Appl. Probab.*, 17(1):102–155, 2007.

[7] J. N. Darroch and E. Seneta. On quasi-stationary distributions in absorbing discrete-time finite Markov chains. *J. Appl. Probab.*, 2:88–100, 1965.

[8] J. N. Darroch and E. Seneta. On quasi-stationary distributions in absorbing continuous-time finite Markov chains. *J. Appl. Probab.*, 4:192–196, 1967.

[9] Charles F. Dunkl and Yuan Xu. *Orthogonal polynomials of several variables*, volume 81 of *Encyclopedia of Mathematics and its Applications*. Cambridge University Press, Cambridge, 2001.

[10] Pablo A. Ferrari, Servet Martínez, and Pierre Picco. Existence of non-trivial quasi-stationary distributions in the birth-death chain. *Adv. Appl. Probab.*, 24(4):795–813, 1992.

[11] David C. Flaspohler. Quasi-stationary distributions for absorbing continuous-time denumerable Markov chains. *Ann. Inst. Statist. Math.*, 26:351–356, 1974.

[12] F. R. Gantmacher. *Matrizentheorie*. Hochschulbücher für Mathematik [University Books for Mathematics], 86. VEB Deutscher Verlag der Wissenschaften, Berlin, 1986. With a foreword by D. P. Želobenko, Translated from the Russian by Helmut Boseck, Dietmar Soyka and Klaus Stengert.

[13] Frédéric Gosselin. Asymptotic behavior of absorbing Markov chains conditional on nonabsorption for applications in conservation biology. *Ann. Appl. Probab.*, 11(1):261–284, 2001.
[14] Göran Högnäs. On the quasi-stationary distribution of a stochastic Ricker model. *Stoch. Process. Appl.*, 70(2):243–263, 1997.

[15] S. Karlin and J. McGregor. Linear growth models with many types and multidimensional Hahn polynomials. In *Theory and application of special functions (Proc. Advanced Sem., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1975)*, pages 261–288. Math. Res. Center, Univ. Wisconsin, Publ. No. 35. Academic Press, New York, 1975.

[16] Harry Kesten. A ratio limit theorem for (sub) Markov chains on $\{1, 2, \dots\}$ with bounded jumps. *Adv. in Appl. Probab.*, 27(3):652–691, 1995.

[17] Kshitij Khare and Hua Zhou. Rates of convergence of some multivariate Markov chains with polynomial eigenfunctions. *Ann. Appl. Probab.*, 19(2):737–777, 2009.

[18] Masaaki Kijima and E. Seneta. Some results for quasi-stationary distributions of birth-death processes. *J. Appl. Probab.*, 28(3):503–511, 1991.

[19] M. Kimura. *The Neutral Theory of Molecular Evolution*. Cambridge University Press, 1983.

[20] Motoo Kimura. Solution of a process of random genetic drift with a continuous model. *Proc. Nat. Acad. Sci.*, 41:144–150, 1955.

[21] J. D. Murray. *Mathematical Biology*. Springer-Verlag, Berlin, second edition, 1993.

[22] Ingemar Nåsell. On the quasi-stationary distribution of the stochastic logistic epidemic. *Math. Biosci.*, 156(1-2):21–40, 1999. Epidemiology, cellular automata, and evolution (Sofia, 1997).

[23] P. K. Pollett. Quasi-stationary distributions: A bibliography. available on http://www.maths.uq.edu.au/~pkp/papers/qsds/, 2011.

[24] E. Seneta and D. Vere-Jones. On quasi-stationary distributions in discrete-time Markov chains with a denumerable infinity of states. *J. Appl. Probab.*, 3:403–434, 1966.
[25] Erik A. van Doorn. Quasi-stationary distributions and convergence to quasi-stationarity of birth-death processes. *Adv. Appl. Probab.*, 23(4):683–700, 1991.

[26] Anton Zettl. *Sturm-Liouville theory*, volume 121 of *Mathematical Surveys and Monographs*. American Mathematical Society, Providence, RI, 2005.