Quasi-stationary distributions and population processes

This survey concerns the study of quasi-stationary distributions with a specific focus on models derived from ecology and population dynamics. We are concerned with the long time behavior of different stochastic population size processes when 0 is an absorbing point almost surely attained by the process. The hitting time of this point, namely the extinction time, can be large compared to the physical time and the population size can fluctuate for large amount of time before extinction actually occurs. This phenomenon can be understood by the study of quasi-limiting distributions. In this paper, general results on quasi-stationarity are given and examples developed in detail. One shows in particular how this notion is related to the spectral properties of the semi-group of the process killed at 0. Then we study different stochastic population models including nonlinear terms modeling the regulation of the population. These models will take values in countable sets (as birth and death processes) or in continuous spaces (as logistic Feller diffusion processes or stochastic Lotka-Volterra processes). In all these situations we study in detail the quasi-stationarity properties. We also develop an algorithm based on Fleming-Viot particle systems and show a lot of numerical pictures.


Introduction
We are interested in the long time behavior of isolated biological populations with a regulated (density-dependent) reproduction. Competition for limited resources impedes these natural populations without immigration to grow indefinitely and leads them to become extinct. When the population's size attains zero, nothing happens anymore and this population's size process stays at zero. This point 0 is thus an absorbing point for the process. Nevertheless, the time of extinction can be large compared to the individual time scale and it is common that population sizes fluctuate for large amount of time before extinction actually occurs. For example, it has been observed in populations of endangered species, as the Arizona ridge-nose rattlesnakes studied in Renault-Ferrière-Porter [52], that the statistics of some biological traits seem to stabilize. Another stabilization phenomenon is given by the mortality plateau. While demographers thought for a long time that the rate of mortality of individuals grows as an exponential function of the age, it has been observed more recently that the rate of mortality slows at advanced ages, or even stabilizes. To capture these phenomena, we will study the long time behavior of the process conditioned on non extinction and the related notion of quasi-stationarity. In particular, we will see that a Markov process with extinction which possesses a quasistationary distribution has a mortality plateau.
In all the following, the population's size process (Z t , t ≥ 0) will be a Markov process going almost surely to extinction. We are interested in looking for characteristics of this process giving more information on its long time behavior. One way to approach this problem is to study the "quasi-limiting distribution" (QLD) of the process (if it exists), that is the limit, as t → +∞, of the distribution of Z t conditioned on non-absorption up to time t. This distribution, which is also called Yaglom limit, provides particularly useful information if the time scale of absorption is substantially larger than the one of the quasi-limiting distribution. In that case, the process relaxes to the quasi-limiting regime after a relatively short time, and then, after a much longer period, absorption will eventually occur. Thus the quasi-limiting distribution bridges the gap between the known behavior (extinction) and the unknown time-dependent behavior of the process. There is another point of view concerning quasi-stationarity. A quasi-stationary distribution for the process (Z t , t ≥ 0) denotes any proper initial distribution on the non-absorbing states such that the distribution of Z t conditioned on non-extinction up to time t is independent of t, t ≥ 0.
There is a large literature on quasi-stationary distributions and Yaglom limits (see for example the large bibliography updated by Pollett [50]) and a lot of references will be given during the exposition. The present paper is by no means exhaustive, but is a survey presenting a collection of tools for the study of QSD concerning specific population's size models. More than the originality of the proofs, we emphasize some general patterns for qualitative and quantitative results on QSD. We also provide a lot of numerical illustrations of the different notions.
In Section 2 of this survey, we will introduce the different notions of QSD and review theoretical properties on QSD and QLD. We will also highlight the relations between QSDs and mortality plateaus. In Section 3, we will study the simple case of QSD for processes in continuous time with finite state space. We develop a simple example to make things more concrete. Thus we will concentrate on QSD for several stochastic population models corresponding to different space and time scales. We will underline the importance of spectral theory as mathematical tool for the research of QSD, in these different contexts. In Section 4, we will consider birth and death processes. We will state results giving explicit conditions on the coefficients ensuring the almost sure extinction of the process, and the existence and uniqueness (or not) of a QSD. We will especially focus on the density-dependence case, when the death rate of each individual is proportional to the population's size (called logistic birth and death process). We will show that in that case, the process goes almost surely to extinction, and that there is a unique QSD, coinciding with the unique QLD. In Section 5, the birth and death process is rescaled by a growing initial population and by small individual weights (small biomass). This process is proved to converge, as the initial population's size tends to infinity, to the unique solution of the deterministic logistic equation, whose unique stable equilibrium is given by the carrying capacity. If the individual birth and death rates are proportional to the population's size while preserving the ecological balance, more numerous are the individuals, smaller are their weights and faster are their birth and death events, reflecting allometric demographies. In that case, the rescaled birth and death process converges, as the initial population size increases, to the solution of a stochastic differential equation with a 1/2-Hölder diffusion coefficient and a quadratic drift, called logistic Feller equation. The existence of the QSD is proved in this case and uniqueness is characterized by a condition meaning the return of the process from infinity in finite time. The proof relies QSD investigation to spectral theory tools developed in a functional L 2 −space. The logistic Feller equation describes the size of a mono-type population where individuals are indistinguable. Motivated by ecological and biodiversity problems, we generalize the model to multi-type populations with intra-specific and inter-specific competition. That leads us to consider stochastic Lotka-Volterra processes. We give conditions ensuring mono-type transient states or coexistence, preserving one or several dominant types in a longer scale. A large place in the survey is given to illustrations obtained by simulations. Some of them derive from an approximation method based on Fleming-Viot interacting particle systems, which is carefully described in Section 6.

A brief bibliography on quasi-stationary distributions
The study of quasi-stationary distributions began with the work of Yaglom on sub-critical Galton-Watson processes [67]. Since then, the existence, uniqueness and other properties of quasi-stationary distributions for various processes have been studied.
In particular, the case of Markov processes on finite state spaces has been studied by Darroch and Seneta, who proved under some irreducibility conditions the existence and uniqueness of the QSD, for both discrete [18] and continuous time settings [19] (detailed proofs and results are reproduced in Section 3 of the present paper). We also refer the reader to the works of van Doorn and Pollett [61] for a relaxation of the irreducibility condition.
The case of discrete time birth and death processes has been treated by Seneta and Vere-Jones [54] and Ferrari, Martínez and Picco [24]. For continuous time birth and death processes, we refer to van Doorn [59]. This last case is quite enlightening, since it leads to examples of processes with no QSD, of processes with a Yaglom limit and an infinite number of QSD and to processes with a Yaglom limit which is the unique QSD (detailed proofs and results are also developed in Section 4 of this survey). For further developments, we may refer to Pakes and Pollett [48] (where results on continuoustime birth and death processes with catastrophic events are obtained), to Bansaye [5] (where a discrete time branching process in random environment is studied), to Coolen-Schrijner [17] (where general discrete time processes are studied) and references therein.
Diffusion processes have also been extensively studied in the past decades, beginning with the seminal work of Mandl [45] for the one-dimensional case and of Pinsky [49] and Gong, Qian and Zhao [29] in the multi-dimensional situation. Martínez, Picco and San Martín [46] and Lladser and San Martín [44] highlighted cases of diffusions with infinitely many quasi-limiting distributions, with a non-trivial dependence on the initial distribution of the process. For recent development of the theory of QSDs for diffusion processes, we refer to Steinsaltz and Evans [57] and Kolb and Steinsaltz [41] where the case of one dimensional diffusions with different boundary conditions is studied. We also emphasize that in the case of Wright-Fisher diffusions and some of its relatives, Huillet [35] derived explicit values of QSDs. Other diffusion processes related to demographic models have been studied in Cattiaux, Collet, Lambert, Martínez, Méléard and San Martín [13], where the case of the Feller logistic diffusion is developed (proofs and results are also written in detail in Section 5 of this paper), and in Cattiaux and Méléard [14], where the case of a two dimensional stochastic Lotka-Volterra system is developed (k types stochastic Lotka-Volterra systems are also studied in Section 5.4 of this survey).
Let us mention the original renewal approach of Ferrari, Kesten, Martínez and Picco [22], also studied recently by Barbour and Pollett [6] in order to provide an approximation method for the QSD of discrete state space Markov processes in continuous time. Other approximation methods have been proposed by Pollett and Stewart [51] and by Hart and Pollett [34]. In this survey, we describe the approximation method based on Fleming-Viot type interacting particle systems, introduced by Burdzy, Holyst, Ingerman et March [10] in 1996 and studied later by Burdzy, Holyst and March [11], Grigorescu and Kang [32], Ferrari and Maric [23], Villemonais [65] [66] and Asselah, Ferrari and Groisman [3].
For studies on the so-called Q-process, which is the process distributed as the original process conditioned to never extinct, we refer the reader to the above cited articles [45], [49], [29], [18], [19] and, for further developments, to the works of Collet, Martínez and San Martín [16] and of Lambert [43] and references therein.

The framework
Let us now introduce our framework in more details. The population's size (Z t , t ≥ 0) is a Markov process taking values in a subset E of N or R + , in a discrete or continuous time setting. If the population is isolated, namely without immigration, then the state 0, which describes the extinction of the population, is a trap. Indeed, if there are no more individuals, no reproduction can occur and the population disappears. Thus if the system reaches 0, it stays there forever, that is, if Z t = 0 for some t, then Z s = 0 for any s ≥ t.
We denote by T 0 the extinction time, i.e. the stopping time We will consider cases for which the process goes almost surely to zero, whatever the initial state is, namely, for all z ∈ E, Before extinction, the process takes its values in the space E * = E\{0}. Any long time distribution of the process conditioned on non-extinction will be supported by E * .
Notations For any probability measure µ on E * , we denote by P µ (resp. E µ ) the probability (resp. the expectation) associated with the process Z initially distributed with respect to µ. For any x ∈ E * , we set P x = P δx and E x = E δx . We denote by (P t ) t≥0 the semi-group of the process Z killed at 0. More precisely, for any z > 0 and f measurable and bounded on E * , one defines For any finite measure µ and any bounded measurable function f , we set and we also define the finite measure µP t by

Definitions, general properties and first examples
There are several natural questions associated with this situation.
Question 1 What is the distribution of the size of a non-extinct population at a large time t ? The mathematical quantity of interest is thus the conditional distribution of Z t defined, for any Borel subset A ⊂ E * , by where ν is the initial distribution of the population's size Z 0 . We will study the asymptotic behavior of this conditional probability when t tends to infinity. The first definition that we introduce concerns the existence of a limiting conditional distribution.
Definition 1. Let α be a probability measure on E * . We say that α is a quasi-limiting distribution (QLD) for Z, if there exists a probability measure ν on E * such that, for any measurable set A ⊂ E * , In some cases the long time behavior of the conditioned distribution can be proved to be initial state independent. This leads to the following definition.
Definition 2. We say that Z has a Yaglom limit if there exists a probability measure α on E * such that, for any x ∈ E * and any measurable set A ⊂ E * , When it exists, the Yaglom limit is a QLD. The reverse isn't true in general and (5) will actually not imply the same property for any initial distribution.
Question 2 As in the ergodic case, we can ask if this Yaglom limit has the conditional stationarity property given by the following definition.
Definition 3. Let α be a probability measure on E * . We say that α is a quasi-stationary distribution (QSD) if, for all t ≥ 0 and any measurable set A ⊂ E * , The main questions are: Does a QSD exists? Is there a unique QSD for the process? We will study examples where QSDs do not exist, or with an infinity of QSDs, or with a unique QSD. The relation between the existence of QSD, QLD and Yaglom limit is clarified in Proposition 1 below. Namely, we will prove that Yaglom limit ⇒ QSD ⇔ QLD.
Question 3 Since the processes we are interested in become extinct in finite time almost surely, the event t < T 0 becomes a rare event when t becomes large. An important question is then to know whether the convergence to the Yaglom limit happens before the typical time of extinction, or if it happens only after very large time periods, in which case the populations whose size are distributed with respect to the Yaglom limit are very rare. Both situations can appear, as illustrated by the simple example of Section 2.3.

Question 4
While most of theoretical results on QLDs, QSDs and Yaglom limits are concerned with existence and uniqueness problems, it would be useful in practice to have qualitative information on the Yaglom limit. We present here particle approximation results and numerical computations of the Yaglom limit for some population's size models, providing some enlightenment on Question 3 above.
Question 5 Another mathematical quantity related to this conditioning is based on a pathwise point of view. In the finite state space case of Section 3 and the logistic Feller diffusion case of Section 5, we will describe the distribution of the trajectories who never attain the trap. This will allow us to define a process, commonly referred to as the Q process for Z. We will prove that the new process defined by this distribution is ergodic, and that its stationary distribution is absolutely continuous with respect to the QSD (but not equal).
The present section is organized as follows. In Subsection 2.1, we state general properties of QLDs, QSDs and Yaglom limits. In Subsection 2.2, we develop the case of the Galton-Watson process. This discrete time process is of historical importance, since the notion of Yaglom limit has originally been developed for this process by Yaglom itself (see [67]). In Subsection 2.3, we develop a very simple example of a process evolving in a finite subset of N. For this process, one can easily prove the existence of the Yaglom limit, the uniqueness of the QSD, and compare the speed of extinction to the speed of convergence to the Yaglom limit. We also provide numerical computation of the relevant quantities.

General properties
Most of the following results are already known by the QSD community. In this section, we emphasize their generality.

QSD, QLD and Yaglom limit
It is clear that any Yaglom limit and any QSD is also a QLD. The reverse implication has been proved by Vere-Jones [63] for continuous time Markov chains evolving in a countable state space. The following proposition extends this result to the general setting.
Proposition 1. Let α be a probability measure on E * . The distribution α is a QLD for Z if and only if it is a QSD for Z.

Remark 1.
When it exists, the Yaglom limit is uniquely defined, while there are processes with an infinity of QSDs (see the birth and death process case of Section 4). We immediately deduce that there exist QSDs which aren't a Yaglom limit. Proof.
(1) If α is a QSD then it is a QLD for Z starting with distribution α.
(2) Assume now that α is a QLD for Z and for an initial probability measure µ on E * . Thus, for any measurable and bounded function f on E * , .
Applying the latter with f (z) = P z (T 0 > s), we get by the Markov property .
Let us now consider f (z) = P z (Z s ∈ A, T 0 > s), with A ⊂ E * . Applying the Markov property again, it yields .
By definition of the QLD α, Pµ(Z t+s ∈A;T 0 >t+s) converges to α(A) and Pµ(T 0 >t+s) Pµ(T 0 >t) converges to P α (T 0 > s), when t goes to infinity. We deduce that, for any Borel set A of E * and any s > 0, α(A) = P α (Z s ∈ A|T 0 > s).
The probability measure α is then a QSD.

Exponential extinction rate
Proposition 2. Let us consider a Markov process Z with absorbing point 0 satisfying (2). Assume that α is a QSD for the process. Then there exists a positive real number θ(α) depending on the QSD such that This theorem shows us that starting from a QSD, the extinction time has an exponential distribution with parameter θ(α) independent of t > 0, given by Proof. By the Markov property, since T 0 ≤ t implies Z t = 0, and P 0 (T 0 > s) = 0. By definition of a QSD, we get Hence we obtain that for all s,t > 0, P α (T 0 > t + s) = P α (T 0 > s)P α (T 0 > t). Let us denote g(t) = P α (T 0 > t). We have g(0) = 1 and, because of (2), g(t) tends to 0 as t tends to infinity. An elementary proof allows us to conclude that there exists a real number θ(α) > 0 such that P α (T 0 > t) = e −θ(α)t .

QSD and exponential moments
The following statement gives a necessary condition for the existence of QSDs in terms of existence of exponential moments of the hitting time T 0 .
Proposition 3 suggests that if the population can escape extinction for too long times with positive probability, then the process has no QSD. This is the case for the critical Galton-Watson process: its extinction time is finite almost surely, but its expectation isn't finite.
Proof. We compute the exponential moment in continuous and discrete time settings. In both cases, it is finite if and only if θ(α) > γ.
In the continuous time setting, (6) says that, under P α , T 0 has an exponential distribution with parameter θ(α). We deduce that, for any θ(α) > γ, In the discrete time setting, (6) says that under P α , T 0 has a geometric distribution with parameter e −θ(α) . We deduce that Since E α (e γT 0 ) is equal to E * E z (e γT 0 )α(dz), the finiteness of the integral implies the last assertion.
Remark 2. In the particular case of an irreducible continuous time Markov chain with state space N such that lim z→+∞ P z (T 0 ≤ t) = 0, ∀t ≥ 0, Ferrari, Kesten, Martínez and Picco [22] proved that the existence of the moment (7) for some z ∈ N and some γ > 0 is equivalent to the the existence of a quasi-stationary distribution.
It is actually not true in any case, as shown by the following counter-example. Let Z be a continuous time random walk on N reflected on 1 and killed at rate 1. Thus, for any λ ∈ [0,1[ and any probability measure µ on N, E µ (e λT 0 ) is finite. Nevertheless the conditional distribution P z (Z t ∈ ·|t < T 0 ) is the distribution of a standard continuous time random walk reflected on 1, which converges to 0 as t tends to infinity. In particular Z has no QLD and thus no QSD.

A spectral point of view
In this section, the results are stated in the continuous time setting. The operator L with domain D(L) denotes the infinitesimal generator of the sub-Markovian semi-group (P t ) associated with the killed process Z. The next proposition links the existence of QSDs for Z and the spectral properties of the dual of the operator L. It is one of the main tools used in a large literature studying QSDs.
Proposition 4. Let α be a probability measure on E * . We assume that there exists a set D ⊂ D(L) such that, for any measurable subset A ⊂ E * , there exists a uniformly bounded sequence (f n ) in D converging point-wisely to 1 A .
Then α is a quasi-stationary distribution if and only if there exists θ(α) > 0 such that We emphasize that the existence of D is always true if the state space E * is discrete. It is also fulfilled if E * is an open subset of R d and if Z is a diffusion with locally bounded coefficients.
Proof. (1) Let α be a QSD for Z. By definition of a QSD, we have, for every Borel set A ⊆ E * , By Theorem 6, there exists θ(α) > 0 such that for each t > 0, We deduce that, for any measurable set A ⊆ E * , αP t (1 A ) = e −θ(α)t α(A), which is equivalent to αP t = e −tθ(α) α. By Kolmogorov's forward equation and by assumption on D, we have In particular, one can differentiate αP t f = E * P t f (x)α(dx) under the integral sign, which implies that α(Lf ) = −θ(α)α(f ), ∀f ∈ D.
By assumption, there exists, for any measurable subset A ⊂ E * , a uniformly bounded sequence (f n ) in D which converges point-wisely to 1 A . Finally, we deduce by dominated convergence that This implies immediately that α is a quasi-stationary distribution for Z.

Long time limit of the extinction rate
Another quantity of interest in the demography and population's dynamics is given by the long time behavior of the killing or extinction rate. In the demography setting, the process Z models the vitality of some individual and t its physical age. Thus T 0 is the death time of this individual. The long time behavior of the extinction rate has been studied in detail by Steinsaltz-Evans [56] for specific cases.
The definition of the extinction rate depends on the time setting: • In the discrete time setting, the extinction rate of Z starting from µ at time t ≥ 0 is defined by r µ (t) = P µ (T 0 = t + 1|T 0 > t).
• In the continuous time setting, the extinction rate of Z starting from µ at time t ≥ 0 is defined by when the derivative exists and is integrable with respect to µ.
Historically (cf. [28]), demographers applied the Gompertz law meaning that this extinction rate was exponentially increasing with time. However in 1932, Greenwood and Irwin [31] observed that in some cases, this behavior was not true. In particular there exist cases where the extinction rate converges to a constant when time increases, leading to the notion of mortality plateau. This behavior of the extinction rate has been observed in experimental situations (see for instance [12]).
The QSDs play a main role in this framework. Indeed, by Proposition 2, if α is a QSD, then the extinction rate r α (t) is constant and given by r α (t) = 1 − e −θ(α) in the discrete time setting θ(α) in the continuous time setting , ∀t ≥ 0.
We refer to the introduction of Steinsaltz-Evans [56] for a nice discussion of the notion of QSD in relationship with mortality plateaus.
In the next proposition, we prove that the existence of a QLD for Z started from µ implies the existence of a long term mortality plateau.
Proposition 5. Let α be a QLD for Z, initially distributed with respect to a probability measure µ on E * . In the continuous time setting, we assume moreover that there exists h > 0 such that L(P h 1 E * ) is well defined and bounded. In both time settings, the rate of extinction converges in the long term: Proof. In the discrete time setting, by the semi-group property and the definition of a QLD, we have The limit is by definition the extinction rate at time 0 of Z starting from α. In the continuous time setting, by the Kolmogorov's backward equation, we have Since L(P h 1 E * ) is assumed to be bounded, we deduce that by the definition of a QLD and by Proposition 4. We also have Finally, we get which allows us to conclude the proof of Proposition 5.

A historical example in discrete time: the Galton-Watson process
The Galton-Watson process is a population's dynamics model in discrete time, whose size (Z n ) n≥0 evolves according to the recurrence formula Z 0 and where (ξ (n) i ) i,n is a family of independent random variables, identically distributed following the probability measure µ on N with generating function g. As defined, Z n is the size of the n th generation of a population where each individual has a random number of children, chosen following µ and independently of the rest of the population. This process has been introduced by Galton and Watson (see [26]) in order to study the extinction of aristocratic surnames. We will assume in the whole section that 0 < µ({0}) + µ({1}) < 1. We denote by m = E(ξ (0) 1 ) the average number of children by individual in our Galton-Watson process. The independence of descendants implies that starting from Z 0 , the process Z is equal to the sum of Z 0 independent Galton-Watson processes issued from a single individual. By this branching property, the probability of extinction for the population starting from one individual is obtained as follows: There are three different situations (see for instance Athreya-Ney [4]): -The sub-critical case m < 1: the process becomes extinct in finite time almost surely and the average extinction time E(T 0 ) is finite.
-The critical case m = 1: the process becomes extinct in finite time almost surely, but E(T 0 ) = +∞.
-The super-critical case m > 1: the process is never extinct with a positive probability, and it yields immediately that E(T 0 ) = +∞.
Theorem 6 (Yaglom [67], 1947). Let (Z n ) n≥0 be a Galton-Watson process with the reproduction generating function g. There is no quasi-stationary distribution in the critical and the super-critical case. In the sub-critical case, the Yaglom limit exists and is the unique QSD of Z. Moreover, its generating functionĝ fulfillŝ g(g(s)) = mĝ(s) + 1 − m, ∀s ∈ [0,1].
Proof. The proof is adapted from Athreya-Ney [4] p. 13-14. In the critical or the supercritical case, we have E 1 (T 0 ) = +∞, which implies that E α (T 0 ) = +∞ for all probability measure α on N * . We deduce from Proposition 3 that there is no QSD.
Since g 1 is convex, Γ is non-decreasing. Moreover m < 1 implies that g n (x) ≥ x, so that g n (s) and 1 −ĝ n (s) are non-decreasing in n. In particular, lim n→∞ĝn (s) exists. Let us denote byĝ(s) its limit and by α the corresponding finite measure (whose mass is smaller than one). In order to prove that α is a probability measure on N * , it is sufficient to prove thatĝ(s) → 1 when s goes to 1. We have Taking the limit on each size, where lim n→∞ Γ(g n (0)) = Γ(1) = m, we deduce that which implies Equation (9). Since lim s→1 g 1 (s) = 1 and m < 1, thenĝ(1) = 1. Finally, α is a QLD for Z starting with distribution ν.
One could think a priori that the functionĝ depends on the starting distribution ν. We prove now that it isn't the case, so that there is a unique QLD, and then a unique QSD, which is also the Yaglom limit of the process (indeed, one could choose ν = δ x , x ∈ N * ). Assume that there exist two generating functionsĝ andĥ which fulfill Equation (9). By induction, we have, for all n ≥ 1 and all s ∈ [0,1], g(g n (s)) = m nĝ (s) + m n−1 + · · · + m + 1 (m − 1), h(g n (s)) = m nĥ (s) + m n−1 + · · · + m + 1 (m − 1).

The simple example of an ergodic process with uniform killing in a finite state space
We present a very simple Markov process with extinction whose quasi-stationary distribution, Yaglom limit, speed of extinction and speed of convergence to the Yaglom limit are very easy to obtain. Let (X t ) t≥0 be an exponentially ergodic Markov process which evolves in the state space E * = {1, · · · ,N }, N ≥ 1. By exponentially ergodic, we mean that there exist a probability measure α on E * and two positive constants C,λ > 0 such that, for all z ∈ {1, · · · ,N } and all t ≥ 0, sup There is no possible extinction for (X t ). Let d > 0 be a positive constant and let τ d be an exponential random time of parameter d independent of the process (X t ). We define the process (Z t ) by setting This model is a model for the size of a population which cannot be extinct, except at a catastrophic event which happens with rate d. Thus we have The conditional distribution of Z t is simply given by the distribution of X t : We deduce that the unique QSD is the Yaglom limit α and that for all z ∈ E * and all Thus in this case, the conditional distribution of Z converges exponentially fast to the Yaglom limit α, with rate λ > 0 and the process becomes extinct exponentially fast, with rate d > 0.
Hence the comparison between the speed of convergence to the Yaglom limit and the speed of extinction will impact the observables of the process before extinction: (a) If λ d, then the convergence to the Yaglom limit happens before the typical time of extinction of the population and the quasi-stationary regime will be observable.
(b) If λ d, then the extinction of the population occurs before the quasi-stationary regime is reached. As a consequence, we are very unlikely to observe the Yaglom limit.
(c) If λ ∼ d, the answer is not so immediate and depends on other parameters, as in particular the initial distribution. Example 1. The population size Z is described by a random walk in continuous time evolving in E = {0,1,2, · · · ,N } with transition rates given by The boundedness of the population size models a constraint of fixed resources which acts on the growth of the population. We will see more realistic fixed resources models including logistic death rate in the next sections. One can check that the quasi-stationary probability measure of Z is given by Numerical simulations. We fix N = 100. In that finite case, one can obtain by numerical computation the whole set of eigenvalues and eigenvectors of the infinitesimal generator L (we use here the software SCILAB). Numerical computation using the fact that λ is the spectral gap of the generator of X gives λ = 0.098. For different values d = 0.001, d = 0.500 and d = 0.098, we compute numerically the mathematical quantities of interest: the extinction probability P z (T 0 > t) = e −dt as a function of t (cf. Figure 1 left picture) and the distance sup i∈E * |P z (Z t = i|T 0 > t) − α({i})| between the conditional distribution of Z t and α as a function of − log P z (T 0 > t), which gives the extinction's time scale. (cf. Figure 1 right picture).
We observe that the convergence to the Yaglom limit happens rapidly in the case ( ) λ = 0.098 d = 0.001. Indeed the distance to the Yaglom limit is equal to 0.05, while the survival probability can't be graphically distinguished from 1. On the contrary, we observe that the convergence happens very slowly in the case (2) λ = 0.098 d = 0.500. Indeed, the distance to the Yaglom limit is equal to 0.05 when the survival probability appears to be smaller than e −15 3 × 10 −7 . The case (·) λ = 0.98 = d is an intermediate case, where the distance to the Yaglom limit is equal to 0.05 when the survival probability appears to be equal to e −3 0.05.
3 The finite case, with general killing's rate

The quasi-stationary distributions
The Markov process (Z t ) t≥0 evolves in continuous time in E = {0,1,...,N }, N ≥ 1 and we still assume that 0 is its unique absorbing state. The semi-group (P t ) t≥0 is the sub-Markovian semi-group of the killed process and we still denote by L the associated infinitesimal generator. In this finite state space case, the operators L and P t are matrices, and a probability measure on the finite space E * is a vector of non-negative entries whose sum is equal to 1. The results of this section have been originally proved by Darroch and Seneta ([18] and [19]).
Theorem 7. Assume that Z is an irreducible and aperiodic process before extinction, which means that there exists t 0 > 0 such that the matrix P t 0 has only positive entries (in particular, it implies that P t has positive entries for t > t 0 ). Then the Yaglom limit α exists and is the unique QSD of the process Z t .
Moreover, denoting by θ(α) the extinction rate associated to α (see Proposition 2), there exists a probability measure π on E * such that, for any i,j ∈ E * , The main tool of the proof of Theorem 7 is the Perron-Frobenius Theorem, which gives us a complete description of the spectral properties of P t and L. The main point is that the matrix P 1 has positive entries. For the proof of the Perron-Frobenius Theorem, we refer to Gantmacher [27] or Serre [55].
Theorem 8 (Perron-Frobenius Theorem). Let (P t ) be a submarkovian semi-group on {1, · · · ,N } such that the entries of P t 0 are positive for t 0 > 0. Thus, there exists a unique positive eigenvalue ρ, which is the maximum of the modulus of the eigenvalues, and there exists a unique left-eigenvector α such that α i > 0 and N i=1 α i = 1, and there exists a unique right-eigenvector π such that π i > 0 and N i=1 α i π i = 1, satisfying In addition, since (P t ) is a sub-Markovian semi-group, ρ < 1 and there exists θ > 0 such that ρ = e −θ . Therefore where A is the matrix defined by A ij = π i α j , and χ > θ and ϑ(e −χt ) denotes a matrix such that none of the entries exceeds Ce −χt , for some constant C > 0.
Proof of Theorem 7. Applying Perron-Frobenius Theorem to the submarkovian semi-group (P t ) t≥0 , it is immediate from (11) that there exists θ > 0 and a probability measure α on E * such that, for any i,j ∈ E * , Summing over j ∈ E * , we deduce that It follows that, for any i,j ∈ E * , Thus the Yaglom limit exists and is equal to α. Since E is finite, we have for any initial We deduce that the Yaglom limit α is the unique QLD of Z, and thus it is its unique QSD. By Proposition 2, we have αP 1 (1 E * ) = e −θ(α) . By (10), this quantity is also equal to e −θ , so that θ = θ(α). The end of Theorem 7 is thus a straightforward consequence of (12) and (13) Remark 3. One can deduce from (12) and (13) that there exists a positive constant C L such that sup where the quantity χ − θ(α) is the spectral gap of L, i.e. the distance between the first and the second eigenvalue of L. Thus if the time-scale χ − θ(α) of the convergence to the quasi-limiting distribution is substantially bigger than the time scale of absorption (χ − θ(α) θ(α)), the process will relax to the QSD after a relatively short time, and after a much longer period, extinction will occur. On the contrary, if χ − θ(α) θ(α), then the extinction happens before the process had time to relax to the quasi-limiting distribution. In intermediate cases, where λ − θ(α) ≈ θ(α), the constant C L , which depends on the whole set of eigenvalues and eigenfunctions of L, plays a main role which needs further investigations.
We generalize Example 1 to a more realistic case where the killing's rate can depend on the size of the population. For instance, it can be higher for a small population than for a big one.
Example 2. Let Z be a Markov process which models a population whose individuals reproduce and die independently, with individual birth rate λ > 0 and individual death rate µ = 1. In order to take into account the finiteness of the resources, the process is reflected when it attains a given value N , that we choose here arbitrarily equal to 100. Thus the process Z evolves in the finite state space {0,1, · · · ,100} and its transition rates are given by i → i + 1 with rate λi, for all i ∈ {1,2, · · · ,99}, i → i − 1 with rate µi, for all i ∈ {1,2,3, · · · ,100}.
The infinitesimal generator of Z is given by , · · · ,99}, L 100,99 = 100 and L 100,100 = −100, The process Z clearly fulfills the conditions of Theorem 7. As a consequence, it has a Yaglom limit α, which is its unique QSD. Moreover, the probability measure α is the unique normalized and positive left eigenvector of L. Since L is a finite matrix of size 100 × 100, one can numerically compute the whole set of eigenvectors and eigenvalues of the matrix (L ij ). This will allow to obtain numerically the Yaglom limit α, its associated extinction rate θ(α), and the speed of convergence χ − θ(α). Moreover, for any t ≥ 0, one can compute the value of e tL , which is equal to P t (the semi-group of Z at time t). Hence, we may obtain the numerical value of the conditioned distribution P Z 0 (Z t ∈ .|t < T 0 ), for any initial size Z 0 . Finally, we are also able to compute numerically the distance between α and the conditioned distribution P Z 0 (Z t ∈ .|t < T 0 ), for any value of λ > 0 and Z 0 ∈ {1, · · · ,100}.
(a) In the first case (λ = 0.9), an individual is more likely to die than to reproduce and we observe that the Yaglom limit is concentrated near the absorbing point 0.
The rate of extinction θ(α) is the highest in this case, equal to 0.100. In fact, the process reaches the upper bound 100 very rarely, so that the behavior of the process is very similar to the one of a linear birth and death process with birth and death rates equal to λ and µ respectively. In Section 4, we study such linear birth and death processes. We show that the Yaglom limit (which exists if and only if λ < µ) is given by a geometric law and θ(α) = µ − λ.
(b) In the second case (λ = µ = 1), we observe that α decreases almost linearly from α 1 to α 100 and the upper bound N = 100 plays a crucial role. In fact, letting N tend to +∞, one would observe that for any i ≥ 1, α i decreases to 0. The extinction rate θ(α) which is equal to 0.014 for N = 100 would also go to 0. The counterpart of this phenomenon for the linear birth and death process studied in Section 4 is that the Yaglom limit will not exist when µ = λ. (c) In the third case (λ = 1.1), the Yaglom limit α is concentrated near the upper bound 100, while the extinction rate is θ(α) = 5.84 × 10 −5 . The comparison with the linear birth and death process is no more relevant, since the important factor in this case is the effect of the upper bound N = 100, which models the finiteness of the resources in the environment.
In Figure 3, we study the effect of the initial position and of the value of the parameter λ on the speed of convergence to the Yaglom limit and on the speed of extinction. We choose the positions Z 0 = 1, Z 0 = 10 and Z 0 = 100, and we look at the two different cases λ = 0.9 and λ = 1.1, which correspond to the subcritical case (a) and to the supercritical case (c) respectively. We represent, for each set of values of (λ,Z 0 ), the distance to the Yaglom limit, sup i∈{1,··· ,100} |P Z 0 (Z t = i|t < T 0 ) − α i | as a function of the time, and the same distance as a function of the logarithm of the survival probability − log P Z 0 (t < T 0 ) (i.e. the extinction time scale). By numerical computation, we also obtain that  In the case (a), we have θ(α) = 0.100 χ − θ(α) = 0.102 and we observe that the speed of convergence depends on the initial position in a non-trivial way: while the survival probability is smaller for the process starting from 10 than for the process starting from 100, the convergence to the Yaglom limit in the extinction's time scale happens faster in the case Z 0 = 10. In the case (c), we have θ(α) = 5.84 × 10 −5 χ − θ(α) = 0.103. The speed of convergence to the Yaglom limit in the extinction's time scale depends on the initial position: if (2) Z 0 = 100, then it is almost immediate; if ( ) Z 0 = 10, the distance between the conditional distribution and the Yaglom limit is equal to 0.05 when the survival probability is around e −0.5 0.61; if (·) Z 0 = 1, then this distance is equal to 0.05 when the survival probability is around e −2.4 0.091.

The Q-process
Let us now study the marginals of the process conditioned to never be extinct.
Theorem 9. Assume that we are in the conditions of Theorem 7. For any i 0 , i 1 , · · · , i k ∈ E * , any 0 < s 1 < · · · , s k < t, the limiting value lim t→∞ Let (Y t , t ≥ 0) be the process starting from i 0 ∈ E * and defined by its finite dimensional distributions Then Y is a Markov process with values in E * and transition probabilities given by It is conservative, and has the unique stationary probability measure (α j π j ) j .
Remark that the stationary probability is absolutely continuous with respect to the QSD, but, contrary to intuition, it is not equal to the QSD.
Proof. Let us denote θ(α) by θ for simplicity. Let i 0 , i 1 , · · · , i k ∈ E * and 0 < s 1 < · · · < s k < t. We introduce the filtration F s = σ(Z u , u ≤ s). Then Thus Let us now show that Y is a Markov process. We have By (15) and Theorem 7, we have Moreover let us compute the infinitesimal generatorL of Y from the infinitesimal generator We thus check that Since Lπ = −θπ, then j∈E * π j L ij = −θπ i and thus j∈E * Lij = 0.

QSD for birth and death processes
We are describing here the dynamics of isolated asexual populations, as for example populations of bacteria with cell binary division, in continuous time. Individuals may reproduce or die, and there is only one child per birth. The population size dynamics will be modeled by a birth and death process in continuous time. The individuals may interact, competing (for example) for resources and therefore the individual death's rate will depend on the total size of the population. In a first part, we recall and partially prove some results on the non-explosion of continuous time birth and death processes. We will also recall conditions on the birth and death rates which ensure that the process goes to extinction in finite time almost surely. In a second part, we concentrate on the cases where the process goes almost surely to zero and we study the existence and uniqueness of quasi-stationary distributions.

Birth and death processes
We consider birth and death processes with rates (λ i ) i and (µ i ) i , that is N-valued pure jump Markov processes, whose jumps are +1 or −1, with transitions where λ i and µ i , i ∈ N, are non-negative real numbers.
Knowing that the process is at state i at a certain time, the process will wait for an exponential time of parameter λ i before jumping to i + 1 or independently, will wait for an exponential time of parameter µ i before jumping to i − 1. The total jump rate from state i is thus λ i + µ i . We will assume in what follows that λ 0 = µ 0 = 0. This condition ensures that 0 is an absorbing point, modeling the extinction of the population. Since these processes have a main importance in the modeling of biological processes, we study in detail their existence and extinction properties, and then their QSDs.
The most standard examples are the following ones.
1. The Yule process. For each i ∈ N, λ i = λi for a positive real number λ, and µ i = 0. There are no deaths. It's a fission model.
2. The linear birth and death process, or binary branching process. There exist positive numbers λ and µ such that λ i = λi and µ i = µi. This model holds if individuals reproduce and die independently, with birth rate equal to λ and death rate equal to µ.
3. The logistic birth and death process. We assume that every individual in the population has a constant birth rate λ > 0 and a natural death rate µ > 0. Moreover the individuals compete to share fixed resources, and each individual j = i creates a competition pressure on individual i with rate c > 0. Thus, given that the population's size is i, the individual death rate due to competition is given by c(i − 1) and the total death rate is µ i = µi + ci(i − 1).
In the following, we will assume that λ i > 0 and µ i > 0 for any i ∈ N * .
We denote by (τ n ) n the sequence of the jump times of the process, either births or deaths. Let us first see under which conditions on the birth and death rates the process is well defined for all time t ≥ 0, i.e. τ = lim n τ n = +∞ almost surely. Indeed, if τ = lim n τ n < ∞ with a positive probability, the process would only be defined for t < τ on this event. There would be an accumulation of jumps near τ and the process could increase until infinity in finite time.
Let us give a necessary and sufficient condition ensuring that a birth and death process does not explode in finite time. The result is already stated in Anderson [2], but the following proof is actually far much shorter and easier to follow. µ k+1 · · · µ n λ k λ k+1 · · · λ n + µ 1 · · · µ n λ 1 · · · λ n .
Proof. 1) Let us more generally consider a pure jump Markov process (X t , t ≥ 0) with values in N, and generator (L ij , i, j ∈ N). We set q i = −L ii . Let (τ n ) n be the sequence of jump times of the process and (U n ) n the sequence of inter-times defined by We also set τ ∞ = lim n→∞ τ n ∈ [0, + ∞]. The process does not explode in finite time almost surely (and is well defined for all time t ∈ R + ), if and only if for each i ∈ N P i (τ ∞ < ∞) = 0.
Let us show that this property is equivalent to the fact that the unique non-negative and bounded solution x = (x i ) i∈N of L x = x is the null solution.
For any i, we set h Indeed, the property is true for n = 0 since i =j L ij q i = 1. Moreover, by conditioning with respect to U 1 and using the strong Markov property, we get since the jump times of the U 1 -translated process are the τ n − U 1 , n ∈ N * . We have By independence of U 1 and X U 1 , we deduce from (16) that As Let (x i ) i be a nonnegative solution of Lx = x bounded by 1, then Since h 1+q i ≥ 0 for all i,j ∈ E, we deduce by iteration from (17) and (18) Let us in the other hand define for any j the quantity z j = E j (e −τ∞ ). Using τ ∞ = lim n τ n , and τ n = n k=1 U k , we deduce by monotone convergence that z j = lim n h (n) j . If the process does not explode a.s., then τ ∞ = ∞ a.s., and lim n h It turns out that the unique nonnegative bounded solution of Lx = x is zero.
If the process explodes with positive probability, then there exists i such that P i (τ ∞ < ∞) > 0. Making n tend to infinity in (17), we get z i = j =i L ij 1+q i z j . Since z i > 0, z is a positive and bounded solution of Lz = z.
2) Let us now apply this result to the birth and death process with λ 0 = µ 0 = 0. Then The equation Lx = x is given by x 0 = 0 and for all n ≥ 1 by Thus, if we set ∆ n = x n − x n−1 , we have ∆ 1 = x 1 and for n ≥ 1, ∆ n+1 = ∆ n µn λn + 1 λn x n . Let us remark that, for any n, ∆ n ≥ 0 and the sequence (x n ) n is nondecreasing. If x 1 = 0, the solution is zero. If not, we get by induction Letting r n = 1 λ n + n−1 k=1 µ k+1 · · · µ n λ k λ k+1 · · · λ n + µ 1 · · · µ n λ 1 · · · λ n , we deduce that r n x 1 ≤ ∆ n+1 ≤ r n x n . Then The boundedness of the sequence (x n ) n is thus equivalent to the convergence of the series k r k .
Corollary 11. Let us consider a BD-process with birth rates (λ i ) i . If there exists a constant λ > 0 such that then the process is well defined on R + .
The proof is immediate. It turns out that the linear BD-processes and the logistic processes are well defined on R + .
Let us now recall under which assumption a BD-process goes to extinction almost surely.
Proof. Let us introduce which is the probability to attain 0 in finite time, starting from i. As before T 0 denotes the extinction time and T I the hitting time of any I. The Markov property yields the induction formula To resolve this equation, we firstly assume that the rates λ i , µ i are nonzero until some fixed level I such that λ I = µ I = 0. We set for each i, u In particular, u (I) 1 = U I 1+U I . Hence, either (U I ) I tends to infinity when I → ∞ and any extinction probability u i is equal to 1 or (U I ) I converges to a finite limit U ∞ and for i ≥ 1, Corollary 13.
1. The linear BD-process with rates λi and µi goes almost surely to extinction if and only if λ ≤ µ.
2. The logistic BD-process goes almost surely to extinction.
Proof. 1) If λ ≤ µ, i.e. when the process is sub-critical or critical, we obtain U I ≥ I − 1 for any I ≥ 1. Then (U I ) I goes to infinity when I → ∞ and the process goes to extinction with probability 1. Conversely, if λ > µ, the sequence (U I ) I converges to µ λ−µ , and an easy computation gives u i = (λ/µ) i .
2) Here we have It is easy to check that (19) is satisfied.

Quasi-stationary distributions for birth and death processes
We consider a BD-process (Z t ) with almost sure extinction. A probability measure α on N * is given by a sequence (α j ) j≥1 of non-negative numbers such that j≥1 α j = 1.
Our first result is a necessary and sufficient condition for such a sequence (α j ) j≥1 to be a QSD for Z. Thereafter we will study the set of sequences which fulfill this condition (we refer the reader to van Doorn [59] for more details).
Theorem 14. The sequence (α j ) j≥1 is a QSD if and only if The next result follows immediately.
Proof of Theorem 14. By Proposition 4 and for a QSD α, there exists θ > 0 such that where L is the infinitesimal generator of Z restricted to N * . Taking the j th component of this equation, we get Summing over j ≥ 1, we get after re-indexing We deduce that θ = µ 1 α 1 , which concludes the proof of Theorem 14.
The study of the polynomials (H n ) has been detailed in Van Doorn [59]. In particular it is shown that there exists a non-negative number ξ 1 such that x ≤ ξ 1 ⇐⇒ H n (x) > 0, ∀n ≥ 1.
We can immediately deduce from this property that if ξ 1 = 0, then there is no quasistationary distribution.
To go further, one has to study more carefully the spectral properties of the semi-group (P t ) and the polynomials (H n ) n , as it has been done in [39], [30] and [59]. From these papers, the polynomials (H n ) n are shown to be orthogonal with respect to the spectral measure of (P t ). In addition, it yields a tractable necessary and sufficient condition for the existence of QSD based on the birth and death rates. The series (S) with general term

plays a crucial role. Remark that (S) converges if and only if
In particular, we obtain 1. If ξ 1 = 0, there is no QSD.
3. If (S) diverges and ξ 1 = 0, then there is a continuum of QSD, given by the one parameter family (α j (x)) 0<x≤ξ 1 : Remark 4. 1. Formula (24) and the approximation method described in Section 6 allow us to deduce a simulation algorithm to get ξ 1 .

Cases with more general birth and death processes have also been studied recently.
Let us mention the infinite dimensional state space setting of Collet, Martínez, Méléard and San Martín [15], where each individual has a type in a continuous state space which influences its birth and death rates, and mutation on the type can occur. The authors give sufficient and quite general conditions for the existence and uniqueness of a QSD. We also refer the reader to the recent work of van Doorn [60], where a transition to the state 0 may occur from any state. The author provides sufficient conditions for the existence of QSDs.
Let us now develop some examples.
The linear case. We assume λ i = λ i ; µ i = µ i and λ ≤ µ. In that case, the BD-process is a branching process, where each individual reproduces with rate λ and dies with rate µ. A straightforward computation shows that the series (S) diverges.
Setting f s : k → s k , we get by the Kolmogorov forward equation,

But the branching property of the process implies
Setting m = 2 λ λ+µ , we deduce that for s < 1, In particular, we deduce that the generating function F t : s → E(s Zt |Z t > 0) of Z t conditioned to Z t > 0 converges when t goes to infinity: We deduce that the Yaglom limit of Z does not exist if λ = µ and is given by the geometric distribution with parameter λ µ if λ < µ: An easy computation yields ξ 1 = µ − λ, since by (24), α 1 = ξ 1 µ . But the series (S) diverges so that for λ < µ, ξ 1 > 0 and there is an infinite number of QSD. If λ = µ, ξ 1 = 0 and there is no QSD.
The logistic case. We assume λ i = λi ; µ i = µi + ci(i − 1). Because of the quadratic term, the branching property is lost and we can not compute the Yaglom limit as above. Therefore, we have no other choice than to study the convergence of the series (S).

We have
Hence the series converges. Thus the Yaglom limit exists and is the unique quasistationary distribution.
One can obtain substantial qualitative information by looking closer to the jump rates of the process. For instance, (λ − µ)/c is a key value for the process. Indeed, given a population size i, the expectation of the next step is equal to i(λ−µ−c(i−1)) λi+µi+ci(i−1) . Then the sign of this expectation depends on the position of i − 1 with respect to (λ − µ)/c: • If i ≤ (λ − µ)/c + 1, then the expectation of the next step will be positive.
We deduce that the region around (λ − µ)/c is stable: it plays the role of a typical size for the population and we expect that the mass of the Yaglom limit is concentrated around it. The value (λ − µ)/c is called (by the biologists) the charge capacity of the logistic BDprocess with parameters λ, µ and c. In the next section, we will consider large population processes, which means logistic BD-processes with large charge capacity.
Example 3. We develop now a numerical illustration of the logistic BD-process case. Across the whole example, the value of the charge capacity λ−µ c is fixed, arbitrarily chosen equal to 9.
In order to illustrate the concept of charge capacity, we represent in Figure 4 a random path of a logistic birth and death process with initial size Z 0 = 1 and with parameters λ = 10, µ = 1 and c = 1. We observe that the process remains for long times in a region around the charge capacity. Moreover, we remark that the process remains mainly below the charge capacity; this is because the jumps rate are higher in the upper region, so that it is less stable than the region below the charge capacity.
Let us now compare the Yaglom limits (numerically computed using the approximation method presented in Section 6) of two different logistic BD processes whose charge capacities are equal to 9 (see Figure 5): (a) Z (a) , whose parameters are λ = 10, µ = 1 and c = 1, , whose parameters are λ = 10, µ = 7 and c = 1/3.   We observe that the Yaglom limits of Z (a) and Z (b) are supported by a region which is around the charge capacity. We also remark that the Yaglom limit of the process Z (b) has a more flat shape than the Yaglom limit of Z (a) . This is because the competition parameter of Z (b) is small in comparison with the birth and death parameters, so that the drift toward the charge capacity is small too, both above and below the charge capacity.
We compute now the distance between the conditioned distribution and the Yaglom limit for the two processes Z (a) and Z (b) for different values of the initial state, namely Z 0 = 1, Z 0 = 10 and Z 0 = 100. The numerical results are represented in Figure 6. We observe a strong dependence between the speed of convergence and the initial position of the processes. In the case of Z (a) , it only takes a very short time to the process starting from 100 to reach the charge capacity, because the competition parameter is relatively high and so is the drift downward the charge capacity. On the contrary, in the case of Z (b) , it takes a longer time for the process to come back from 100 to the charge capacity, so that the speed of convergence to the Yaglom limit is slow. In both cases, the convergence to the Yaglom limit happens very fast when starting from the value 10, because it is near the charge capacity.

A large population model
We are now rescaling logistic birth and death processes with a parameter K ∈ N * modeling a large population with small individuals, i.e. a population with a large initial size of order K and a large charge capacity assumption. The individual's weights (or biomasses) are assumed to be equal to 1 K and we study the limiting behavior of the total biomass process ( Z K t K , t ≥ 0) when K tends to infinity, Z K t being the population's size at time t. In what follows, λ, µ and c are fixed positive constants. In Subsection 5.1.1, individual birth and death rates are assumed to be constant and the competition rate depends linearly on the individual biomass 1 K . In Subsection 5.1.2, we investigate the qualitative differences of evolutionary dynamics across populations with allometric demographics: life-lengths and reproduction times are assumed to be proportional to the individual's weights.
In both cases, the charge capacity of (Z K ) will be (λ − µ)K/c.

Convergence to the logistic equation
Given a parameter K scaling the population's size, we consider the logistic BD-process Z K with birth, death and competition parameters λ, µ and c/K respectively. We assume that the initial value of Z K is of order K, in the sense that there exists a non-negative real random variable X 0 such that We consider the total biomass process defined by X K = Z K /K for all K ≥ 1 and are interested in the limit of X K when K → ∞. The transitions of the process (X K t , t ≥ 0) are the following ones: Theorem 17. Assume that X 0 is a positive number x 0 . Then, the process (X K t , t ≥ 0) converges in law in D([0,T ], R + ) to the unique continuous (in time) deterministic function solution of

Remark 5.
The function x is thus solution of the ordinary differential equatioṅ called logistic equation. This equation has been historically introduced as the first macroscopic model describing populations regulated by competition (cf. [64]). In Theorem 17 above, it is obtained as the limit of properly scaled stochastic jump models.
The function x solution of (25) hits 0 in finite time if λ < µ, while it remains positive forever if λ > µ, converging in the long term to its charge capacity λ−µ c . Thus at this scale extinction does not happen.
Proof of Theorem 17. The Markov process (X K t , t ≥ 0) is well defined and its infinitesimal generator is given, for any measurable and bounded function φ, by Hence, by Dynkin's theorem ( [21] Prop. IV-1.7), the process defines a local martingale, and a martingale, as soon as each term in (27) is integrable.
In particular, taking φ(x) = x gives that (X K t , t ≥ 0) is a semimartingale and there exists a local martingale M K such that Since x 0 is deterministic and using a localization argument, we deduce that E(sup t≤T (X K t ) 2 ) < ∞. Moreover, taking φ(x) = x 2 applied to (27), and comparing with Itô's formula applied to (X K ) 2 prove that (M K ) is a square-integrable martingale with quadratic variation process Let us now study the convergence in law of the sequence (X K ), when K tends to infinity. For any K, the law of X K is a probability measure on the trajectory space D T = D([0,T ], R + ), namely the Skorohod space of left-limited and right-continuous functions from [0,T ] into R + , endowed with the Skorohod topology. This topology makes D T a Polish state, that is a metrizable complete and separable space, which is not true if D T is endowed with the uniform topology. See Billingsley [9] for details.
The proof of Theorem 17 is obtained by a compactness-uniqueness argument. The uniqueness of the solution of (25) is immediate.
By a natural coupling, one may bound the birth and death process X K stochastically from above by the Yule process Y K started from x 0 , which jumps from x to x + 1 K , at the same birth times than X K . One easily shows that sup K E(sup t≤T (Y K t ) 3 ) < ∞ and thus From this uniform estimate, we deduce the uniform tightness of the laws of X K (as probability measures on D T ), using the Aldous criterion (cf. Aldous [1], Joffe-Métivier [37]). Then, by Prokhorov's Theorem, the compactness of the laws of (X K ) follows. Since sup t≤T |X K t − X K t − | ≤ 1 K and the function x → sup t≤T |x t − x t − | is continuous on D T , each limiting value (in law) of the sequence (X K ) will be a pathwise continuous process. In addition using (29) and (5.1.1), it can be shown that lim K→∞ E( M K t ) = 0. Then, the random fluctuations disappear when K tends to infinity and the limiting values are deterministic functions. Now it remains to show that these limiting values are solutions of (25), which can be done similarly to the proof of Theorem 18 (4) stated below.

The logistic Feller diffusion process
In this section, we study the logistic BD-processes Z K with birth and death rates given by γ K + λ and γ K + µ respectively. Here γ, λ and µ are still positive constants. We assume that the competition parameter is given by c/K, so that the charge capacity of Z K is still (λ − µ)K/c. Remark 6. This BD-process (Z K t ) t can also be interpreted as a time-rescaled BD-process Y K Kt , whose birth, death and competition parameters are given by γ + λ/K, γ + µ/K and c/K 2 respectively, that is a critical BD-process with small pertubations.
We are considering as in Section 5.1.1 the sequence of processes X K defined for all t ≥ 0 by The transitions of the process (X K ) are given by Formula (28) giving the semi-martingale decomposition of X K will stay true in this case with another square integrable martingale part N K such that One immediately observes that the expectation of this quantity does not tend to zero as K tends to infinity. Hence the fluctuations will not disappear at infinity and the limit will stay random. Let us now state the convergence theorem.
Theorem 18. Assume γ, c, λ, µ > 0 and λ > µ. i) Assume that the sequence (X K 0 ) K converges in law to X 0 with E(X 3 0 ) < ∞. Then the sequence of processes (X K ) K with transitions (30) converges in law in P(D T ) to the continuous process X, defined as the unique solution of the stochastic differential equation where (B t ) t∈[0,+∞[ is a standard Brownian motion.
ii) Let us introduce for each y ≥ 0 the stopping time For any x ≥ 0, we get P x (T 0 < ∞) = 1.
When c = 0, Equation (31) defines the super-critical Feller diffusion process and will explode with positive probability. In the general case where c = 0, it defines the so-called logistic Feller diffusion process following the terminology introduced by Etheridge [20] and Lambert [42]. Let us remark that the quadratic term driven by c regulates the population size, which fluctuates until it attains the absorbing point 0. Theorem 18 shows that an accumulation of a large amount of birth and death events may create stochasticity, often called by biologists ecological drift or demographic stochasticity. Contrarily to the previous case (Theorem 17), the limiting process suffers extinction almost surely.
Proof. As for Theorem 17, the proof is based on a uniqueness-compactness argument.
(1) The uniqueness of the solution of (31) follows from a general existence and pathwise uniqueness result in Ikeda-Watanabe [36] Section IV-3 or Karatzas-Shreve [38]. For a stochastic differential equation with σ and b smooth enough, the existence and pathwise uniqueness are determined thanks to the following scale functions: for x > 0, More precisely, it is proved that (i) ∀x > 0, P x (T 0 < T ∞ ) = 1 and (ii) Λ(+∞) = +∞ ; κ(0 + ) < +∞ are equivalent. In that case, pathwise uniqueness follows, and then uniqueness in law. In our situation, the coefficients are given by so that the functions Λ and κ satisfy (ii). Thus the SDE (31) has a pathwise unique solution which reaches 0 in finite time almost surely.
(2) Let us assume that E(X 3 0 ) < ∞ and prove that sup K E(sup t≤T (X K t ) 2 ) < ∞. The infinitesimal generator of X K is given bỹ With φ(x) = x 3 , we obtain that where M K is a local martingale. Using a standard localization argument and we get where C is independent of K. Gronwall's lemma yields Now, thanks to (35) and Doob's inequality, we deduce from the semi-martingale decomposition of (X K t ) 2 (obtained using φ(x) = x 2 ), that (3) As previously, the uniform tightness of the laws of (X K ) is obtained from (36) and the Aldous criterion [1]. Therefore, the sequence of laws is relatively compact and it remains to characterize its limit values.
(4) As in the proof of Theorem 17, we remark that the limiting values only charge the set of continuous trajectories, since sup t≤T |∆X K t | ≤ 1 K . Let Q ∈ P(C([0,T ], R + )) be a limiting value of the sequence of laws of the processes X K . We will identify Q as the (unique) law of the solution of (31) and the convergence will be proved. Let us denote C T = C([0,T ], R + ) and define, for φ ∈ C 2 b and t > 0, the function which is continuous Q-a.s.. Let us show first that the process (ψ t (X)) t is a Q-martingale.
For x ∈ R + , we define Using Taylor's expansion, we immediately get (withL K defined in (34)) where C doesn't depend on x and K. By (36), we deduce that E |L K φ(X K t ) − Lφ(X K t )| tends to 0 as K tends to infinity, uniformly for t ∈ [0,T ].
For s 1 < · · · < s k < s < t, for g 1 , · · · , g k ∈ C b , we introduce the function H defined on the path space by

Let us show now that
which will imply that (ψ t (X)) t is a Q-martingale.
In another way, this quantity is equal to

The first term is equal to
ds and tends to 0 by (36) and (37). The second term is equal to E(H(X K ) − H(X)). The function X → H(X) is continuous and since H(X) ≤ C 1 + t s (1 + X 2 u )du , it is also uniformly integrable by (35). This leads the second term to tend to 0 as K tends to infinity. Therefore, it turns out that (38) is fulfilled and the process ψ t (X) = φ(X t ) − φ(X 0 ) − t 0 Lφ(X s )ds is a Q-martingale. By (36) and taking φ(x) = x leads to X t = X 0 + M t + t 0 ((λ − µ)X s − cX 2 s )ds, where M is a martingale. Taking φ(x) = x 2 on the one hand and applying Itô's formula for X 2 t on the other hand allow us to identify By the representation theorem proved in [38] Theorem III-4.2 or in [36], there exists a Brownian motion B such that That concludes the proof.

Statement of the results
We are now interested in studying the quasi-stationarity for the logistic Feller diffusion process solution of the equation where the Brownian motion B and the initial state Z 0 are given, and r and c are assumed to be positive. (We have assumed that γ = 1/2). The results and proofs that are presented in Section 5.2 have been obtained by Cattiaux, Collet, Lambert, Martínez, Méléard and San Martín [13].
Let us firstly state the main theorem of this part.
Theorem 19. Assume that Z 0 , r and c are positive. Then the Yaglom limit of the process Z exists and is a QLD for Z starting from any initial distribution. As a consequence, it is the unique QSD of Z.
Remark 7. 1) The theory studying the quasi-stationary distributions for one-dimensional diffusion processes started with Mandl [45] and has been developed by many authors. See in particular [16], [47], [57], [40]. Nevertheless in most of the papers, the diffusion and drift coefficients are regular and the "Mandl's condition" κ(+∞) = ∞ (see (33)) is assumed. This condition is not satisfied in our case because of the degeneracy of the diffusion and the unboundedness of the drift coefficient.
2) Theorem 19 differs from the results obtained in case of drifts going slower to infinity. For example, Lambert [43] proves that if c = 0 and r ≤ 0, then either r = 0 and there is no QSD, or r < 0 and there is an infinite number of QSD. Lladser and San Martín [44] show that in the case of the Ornstein-Uhlenbeck process dY t = dB t − Y t dt, killed at 0, there is also a continuum of QSD. In the logistic Feller diffusion situation as in the logistic BD-process, the uniqueness comes from the quadratic term c X 2 t induced by the ecological constraints.
3) We have seen that the rescaled charge capacity of the logistic birth and death process converges to the charge capacity of the logistic Feller diffusion. However, whether the rescaled Yaglom limit of the logistic birth and death process converges to the Yaglom limit of the logistic Feller diffusion process remains an open problem.
In order to prove Theorem 19, we firstly make a change of variable and introduce the process (X t , t ≥ 0) defined by X t = 2 √ Z t . Of course, X is still absorbed at 0 and QSDs for Z will be easily deduced from QSDs for X. From now on, we focus on the process (X t ).
An elementary computation using Itô's formula shows that (X t ) is the Kolmogorov diffusion process defined by with Mention that the function q is continuous on R * + but explodes at 0 as 1 2x and at infinity as c 8 x 3 . The strong (cubic) downward drift at infinity will force the process to live essentially in compact sets. That will provide the uniqueness of the QSD, as seen below.
We introduce the measure µ, defined by where Q is given by In particular −Q/2 is a potential of the drift −q. The following result clearly implies Theorem 19.
Theorem 20. [13] Assume that X 0 , r and c are positive. Then the Yaglom limit α of the process X exists. Moreover, there exists a positive function η 1 ∈ L 2 (dµ) such that 1.
4. the QSD α attracts all initial distribution, which means that α is a QLD for X starting from any initial distribution.
The proof of Theorem 20 will be decomposed in the next subsections.

Spectral theory for the killed semi-group
As previously we are interested in the semi-group of the killed process, that is, for any x > 0, t > 0 and any with the associated infinitesimal generator given for φ ∈ C 2 c ((0, + ∞)) by We are led to develop a spectral theory for this generator in L 2 (µ). Though the unity function 1 does not belong to L 2 (µ), this space is the good functional space in which to work. The key point we firstly show is that, starting from x > 0, the law of the killed process at time t is absolutely continuous with respect to µ with a density belonging to L 2 (µ). The first step of the proof is a Girsanov Theorem.

Proposition 21. For any bounded Borel function
where E Wx denotes the expectation with respect to the Wiener measure starting from x and ω the current point in Ω.
Proof. It is enough to show the result for non-negative and bounded functions F . Let ε ∈ (0,1) and τ ε = T ε ∧ T 1/ε . Let us choose some ψ ε which is a non-negative C ∞ function with compact support included in ]ε/2, 2/ε[ such that ψ ε (u) = 1 if ε ≤ u ≤ 1/ε. For all x such that ε ≤ x ≤ 1/ε the law of the diffusion (39) coincides up to τ ε with the law of a similar diffusion process X ε obtained by replacing q with the cutoff function q ε = qψ ε . For the latter we may apply Novikov criterion (cf. [53] p.332), ensuring that the law of X ε is given via Girsanov's formula. Hence integrating by parts the stochastic integral. But 1 t<τε is non-decreasing in ε and converges almost surely to 1 t<T 0 both for W x and for P x (since P x (T 0 < ∞) = 1)). Indeed, almost surely, lim ε→0 X τε = lim ε→0 X τε = lim ε→0 ε = 0 so that lim ε→0 τ ε ≥ T 0 . But τ ε ≤ T 0 yielding the equality. It remains to use Lebesgue monotone convergence theorem to finish the proof.

Theorem 22.
For all x > 0 and all t > 0 there exists a density function r(t,x,.) that satisfies for all bounded Borel function f . In addition, for all t > 0 and all x > 0, Proof. Define Denote by e −v(t,x,y) = (2πt) − 1 2 exp − (x−y) 2 2t the density at time t of the Brownian motion starting from x. According to Proposition 21, we have In other words, the law of X t restricted to non extinction has a density with respect to µ given by r(t,x,y) = E Wx (G|ω t = y) e −v(t,x,y)+Q(y) .
Thanks to Theorem 22, we can show, using the theory of Dirichlet forms (cf. Fukushima's book [25]) that the infinitesimal generator L of X, defined by (5.2.2), can be extended to the generator of a continuous symmetric semi-group of contractions of L 2 (µ) denoted by (P t ) t≥0 . In all what follows, and for f, g ∈ L 2 (µ), we will denote f, g µ = R + f (x)g(x)µ(dx). The symmetry of P t means that P t f, g µ = f, P t g µ . In Cattiaux et al. [13], the following spectral theorem in L 2 (µ) is proved.
The proof of this theorem is based on a relation between the Fokker-Planck operator L and a Schrödinger operator. Indeed, let us set for g ∈ L 2 (dx), P t g = e −Q/2 P t (g e Q/2 ). P t is a strongly semi-group on L 2 (dx) with generator defined for g ∈ C ∞ c ((0, + ∞)) bỹ The spectral theory for such Schrödinger operator with potential (q 2 −q ) 2 on the line (or the half-line) is well known (see for example the book of Berezin-Shubin [7]), but the potential (q 2 −q ) 2 does not belong to L ∞ loc as generally assumed. Nevertheless, in our case inf(q 2 − q ) > −∞, which ensures the compactness of the operatorsL andP t .
The following corollary of Theorem 23 is a generalization of the Perron-Frobenius Theorem in this infinite-dimensional framework.
Corollary 24. For any bounded and measurable function f , we have Proof. Fix t > 0 and let f be a bounded measurable function on R * + . Let us first prove that P t f belongs to L 2 (µ). On the one hand, we have On the other hand, by Proposition 21, we have, for all x ∈ R * + , But the function we deduce that there exists a constant K t > 0 independent of x and f such that Finally (P t f ) 2 is integrable with respect to µ, so that P t f ∈ L 2 (µ).

Now we deduce from Theorem 23 that
If f belongs to L 2 (µ), then the symmetry of P t implies that Since η i ∈ L 1 (µ), we deduce from the Dominated Convergence Theorem that the equality P t f, η i µ = e −λ i t f, η i µ extends to all measurable bounded functions. This and the equality (45) allow us to conclude the proof of Corollary 24.
In particular, h : y → E * r(1,x,y)ν(dx) belongs to L 2 . Then we deduce from (46) that We conclude that α attracts any compactly supported probability measure.
Let us now prove that α attracts all initial distributions ν supported in (0,∞). We want to show that, for any probability measure ν on R * + , for any Borel set A, we get This is part (4) of Theorem 20 and it clearly implies the uniqueness of the QSD for X (and hence for Z).
Proposition 25. For any a > 0, there exists y a > 0 such that sup x>ya E x (e aTy a ) < ∞.
Proving that α attracts all initial distributions requires the following estimates near 0 and ∞.
For the second limit, we set A 0 := sup where y λ 1 is taken from Proposition 25. Then for large M > y λ 1 , we have Using lim u→∞ e λ 1 u P x 0 (T 0 > u) = η 1 (x 0 ) η 1 ,1 µ , we obtain B 0 := sup u≥0 e λ 1 u P x 0 (T 0 > u) < ∞. Then Let ν be any fixed probability distribution whose support is contained in (0,∞). We must prove (47). We begin by claiming that ν can be assumed to have a strictly positive density h, with respect to µ. Indeed, let which implies that r(1,x,y)ν(dx) is finite dy−a.s.. Finally, define h = / dµ. Notice that for dρ = hdµ showing the claim.
Consider M > ε > 0 and any Borel set A included in (0,∞). Then is bounded by the sum of the following two terms .
We have the bound Thus, from Lemma 26 we get On the other hand we have since α attracts any compactly supported probability measures , and the result follows.
Example 4. We develop now a numerical illustration of this logistic Feller diffusion case. As for the logistic birth and death process (see Example 3, Section 4), the value of the charge capacity r c will remain equal to the fixed value 9 across the whole example. We begin by showing in Figure 7 a random path of a logistic Feller diffusion process with initial size Z 0 = 1 and with parameters r = 9 and c = 1 (an Euler method is used for the numerical simulation of the random path). We observe that the process quickly attains the value of the charge capacity and remains around it for a long time.
We compare now the Yaglom limits of two different logistic Feller diffusion processes whose charge capacity is equal to 9 (see Figure 8):  (a) Z (a) , whose parameters are r = 9 and c = 1, , whose parameters are r = 3 and c = 1/3.
As for the logistic BD-processes, we observe that the two Yaglom limits are centered around the charge capacity. But as a consequence of the relatively weak noise around the charge capacity, the Yaglom limit has clearly a smaller variation around this value in the logistic Feller diffusion case than in the logistic BD process case. We also observe that the smaller are the parameters, the flatter is the Yaglom limit and with a similar explanation as in the logistic BD-process case.
We observe now the distance between the conditional distributions of Z (a) and Z (b) and their respective Yaglom limits, for different initial states, namely Z 0 = 1, Z 0 = 10 and Z 0 = 100. The results, computed with the help of the approximation method studied in Section 6, are represented on figure 9. For both Z (a) and Z (b) , the speed of convergence to the Yaglom limit is the highest for Z 0 = 10, which is quite intuitive since the value of the charge capacity is 9. We also observe that it is higher for the processes starting from 100 than for the processes starting from 1. In particular, this behavior is different than in the logistic birth and death process case.

The Q-process
Let us now describe the law of the trajectories conditioned to never attain 0. Theorem 27. [13] Let us fix a time s and consider B a measurable subset of C([0,s],R + ). Then for any x ∈ R * + , where Q x is the law of a continuous process with transition probabilities given by q(s,x,y)dy, with q(s,x,y) = e λ 1 s η 1 (y) η 1 (x) r(s,x,y) e −Q(y) .

The case of a multi-type population
Until now, we have considered a population where all individuals have the same ecological parameters. This biological assumption corresponds to the case where individuals have the same type. In this section, we generalize the previous study to a population composed of k different types. The population size process describing the dynamics of each subpopulation is given by a k-dimensional stochastic Lotka-Volterra process Z = (Z 1 t , · · · , Z k t ) t≥0 (SLVP), which describes the size of a k-types density dependent population. This model generalizes to k types the 2-types density dependent model introduced by Cattiaux and Méléard [14].
More precisely, we consider for i,j ∈ {1, · · · ,k} the coefficients The process Z takes its values in (R + ) k and is solution of the stochastic differential system where (B i ) i=1,··· ,k are independent standard Brownian motions independent of the initial data Z 0 . The system (50) can be obtained as (31) as approximation of renormalized k-types birth and death processes in case of large population and small life lengths and reproduction times. The coefficients r i are the asymptotic growth rates of i-type's populations. The positive coefficients γ i can be interpreted as demographic parameters describing the ecological timescale. The coefficient c ij , for i,j = 1, · · · ,k, represents the pressure felt by an individual holding type i from an individual with type j. Intra-specific competition is modeled by the rates c ii , while inter-specific competition is described by the coefficients c ij > 0, i = j. If c ij = 0 for all i = j, the stochastic k-dimensional process reduces to k independent Feller logistic diffusion processes. Extinction of the population is modeled by the absorbing state (0, · · · , 0) and the extinction of the subpopulation of type i is modeled by the absorbing set We denote by D the open subset of R k defined by D = (R * + ) k and by ∂D its boundary. We denote by T 0 the first hitting time of (0, · · · ,0), by T A the first hitting time of some subset A and thus by T ∂D the exit time of D. Of course, some of these stopping times are comparable. For example if the initial condition belongs to D, On the other hand, T H i and T H j are not directly comparable for i = j.
Let us prove the existence of the SLVP .
Proposition 29. The process (Z t ) t is well defined on R + . In addition, for all x ∈ (R + ) k , P x (T 0 < +∞) = 1 and there exists λ > 0 such that Let us now state the first theorem, which is concerned by conditioning on the co-existence of the k types.
Theorem 32. Under the balance conditions (55), there exists a unique quasi-stationary distribution ν for the process (X) and the absorbing set ∂D, which is the quasi-limiting distribution starting from any initial distribution: for any µ on D and any A ⊂ D, Furthermore, there exist λ > 0 and a positive function η such that Proof. The proof of the existence of a quasi-stationary distribution results from the spectral theory for the semi-group of the killed process (P t ) (related to X) established in Cattiaux-Méléard [14], Appendices A, B, C. Define the reference measure on (R + ) k by As in Subsection 5.2.2, one builds a self-adjoint operator on L 2 (µ) which coincides with P t for bounded functions belonging to L 2 (µ). Its generator L is self-adjoint on L 2 (µ) and We check that the assumptions required in [14] Theorem A.4 are satisfied and therefore, the operator −L is proved to have a purely discrete spectrum of non-negative eigenvalues and the smallest one λ is positive. The corresponding eigenfunction η is proved to be in L 1 (µ) and the probability measure ν = η dµ D η dµ is the Yaglom limit. let us emphasize that the uniqueness of the quasi-stationary distribution results by [14] Proposition B.12 from the ultracontractivity of the semi-group P t (ultracontractivity means that P t maps continuously L 2 (µ) in L ∞ (µ) for any t > 0). The proof of the latter is easily generalized from the two-types case ( [14] Proposition B.14) to the k-types case.
Theorem 32 shows that in some cases, a stabilization of the process with co-existence of the k types will occur before one of these types disappears. Let us now come back to our initial question: the long-time behavior of the process conditioned on non-extinction. For each i = 1, . . . ,k, we denote by λ i the smallest eigenvalue related to the purely discrete spectrum of the generator for the i-axis diffusion defined by the stochastic differential equation (52).
Theorem 33. Under the balance conditions (55), there exists a Yaglom limit m for the process (X) conditioned on non extinction: for any x = 0, for any A ⊂ D, The support of this measure is included in the k axes. Furthermore, if there exist i 1 ,..., i l ∈ {1, · · · , k} such that λ i 1 = · · · = λ i l < min i =i 1 ,...,i l λ i , then this QSD is concentrated on the axes of coordinates i 1 ,..., i l .
Remark that by Proposition 31, the process really leaves ∂D by the interior of H i . Each system (U (i),j ) j =i is a (k−1)-type kolmogorov process (53) with balance conditions. Hence, by our induction assumption (A k−1 ), there exist for each i ∈ {1, · · · , k} a positive constant v i , a positive function η i and a probability measure ν i on H i such that (56) holds for (U (i),j ) j =i , i ∈ {1, · · · ,k}.

Let us define
For any bounded measurable function f on (R + ) k such that f (0) = 0 and for all t ≥ 0, we have where we used the fact that X (k) reaches ∂D by hitting the interior of one and only one H i . By Theorem 32, there exist a positive constant λ , a positive function η and a probability measure ν on (R * + ) k such that Moreover, a similar coupling argument as in [14] yields λ > v min . We deduce that For each i ∈ {1, · · · , k}, we have by the Markov property By the induction assumption is uniformly bounded. Moreover the inequality 0 < v i < λ and Proposition 3 ensure that E x e v min T ∂D < +∞. Using the convergence property of the induction assumption (A k−1 ) and the dominated convergence theorem, we deduce that We have then which gives us the first part of the induction assumption (A k ).
In order to prove the second part of (A k ), let us introduce the SLVP Y (k) with coefficients (c ij ) defined by c kk = c kk , c ij = c ij and c ki = c ik = 0, ∀i,j = 1, · · · ,k − 1.
By the same coupling argument as above, the return time to ∂D for X (k) is stochastically dominated by the return time to ∂D for Y (k) , i.e. P x (X Since the k − 1 first components of Y (k) are independent of the last one and since . On the one hand, the dynamic of (Y (k),1 , · · · ,Y (k),k−1 ) is the same as U (k) , so that, by the second part of the induction assumption (A k−1 ) and by the definition of v min , sup t≥0,x∈D e v min t P x ((Y (k),1 t , · · · ,Y (k),k−1 t ) ∈ (R * + ) k−1 ) < +∞.
On the other hand, Y (k),k is a one dimensional SLVP, thus we deduce from (A 1 ) that there exists a positive constant λ 1 such that sup t≥0,x∈D e λ 1 t P x (Y (k),k t ∈ R * + ) < +∞.
As a consequence, we have sup t≥0,x∈D e (v min +λ 1 )t P x (X (k) t ∈ D) ≤ sup t≥0,x∈D e (v min +λ 1 )t P x (Y (k) t ∈ D) < +∞ and we deduce that sup x∈D E x (e v min T ∂D ) < +∞.
For any bounded measurable function f , this immediately leads us to sup t≥0,x∈E Moreover, by Equality (58) and the second part of (A k−1 ), we deduce that, for each i ∈ {1, · · · ,k}, sup By Equality (57), the second part of the induction assumption (A k ) is thus proved.
By induction on k ≥ 1, we conclude that Assumption (A k ) is true for any k ≥ 1, thus Theorem 33 follows. (c) the probability P (1,1,1) (Z i t > 0, Z j t = 0, Z k t = 0 | T 0 > t) of existence of one and only one type i, for each type i = 1, 2 and 3.
As expected, the 3-type mode disappears quickly and the 2-type modes are transient. We also observe that the probability P (1,1,1) (Z 1 t > 0, Z 2 t = 0, Z 3 t = 0 | T 0 > t) converges to 1 when t increases, meaning that the last state of the population before extinction is monotype with type 1. It turns out that the support of the conditional law P (1,1,1) ((Z 1 t ,Z 2 t ,Z 3 t ) ∈ · | T 0 > t) becomes more and more concentrated on R * + × {0} × {0} in the long time. The Yaglom limit is thus equal to ν 1 ⊗ δ (0,0) , where ν 1 is the Yaglom limit of the process absorbed at 0 and is represented in Figure 11.
6 Simulation: the Fleming-Viot system As seen in the previous sections, the spectral theory is a powerful tool to prove existence and eventually uniqueness of a QSD for a given process Z. It is based on the equivalence property of Proposition 4, stating that a probability measure α on E * is a QSD for the killed process Z if and only if αL = −θ(α)α, where L denotes the infinitesimal generator of Z and θ(α) a positive constant. In some cases, such as in the finite state space case, one can easily compute numerically the whole set of eigenvalues and eigenvectors of L as seen in Example 1 and Example 2. For these numerical illustrations, we used the software SCILAB and its function spec. We also refer to [62] for a detailed description of some algorithms available in MATLAB for the computation of eigenfunctions and eigenvalues in large (but finite) state space cases.
In other cases, such as the logistic birth and death process of Section 4 and the logistic Feller diffusion of Section 5, solving numerically Equation (59) is too hard and we use a different approach. This approach consists in approximating the QSD and the conditioned distribution P z (Z t ∈ .|t < T 0 ) by the empirical distribution of a simulable interacting particle system. This Fleming-Viot type system, built for any number of particles N ≥ 2, has been introduced by Burdzy, Holyst and March [10] and explored in [11] and in Grigorescu-Kang [32] for d-dimensional killed Brownian motions. It has also been studied in Villemonais [65] for multi-dimensional diffusion processes with unbounded drifts and a general result is available in [66]. Similar systems have also been considered by Ferrari-Maric [23] for continuous Markov chains in a countable state space. In this section, we explain the approximation method based on the Fleming-Viot type interacting particle systems.
Let Z be a killed Markov process which evolves in the state space E. Fix N ≥ 2 and let Z 0 ∈ E be its initial value. The interacting particle system with N particles (Z 1 , · · · ,Z N ) starts from (Z 0 , · · · ,Z 0 ) and belongs to (E * ) N . The particles evolve independently from this initial position according to the law of the killed Markov process Z, until one of them hits the state 0. At that time τ 1 , the killed particle jumps to the position at τ 1 of one of the N − 1 remaining particles, chosen uniformly among them. Then the particles evolve independently according to the law of Z until one of them attains 0 (time τ 2 ), and so on. The sequence of jumps is denoted by (τ n ) n and we set τ ∞ = lim n→∞ τ n .
This procedure defines the (E * ) N -valued process (Z 1 , · · · , Z N ) for all time t ∈ [0,τ ∞ [. Figure 12 shows an illustration of such a system with two particles evolving between their jumps as Markov processes absorbed in 0 and 1.
If τ ∞ = +∞ almost surely, then the Fleming-Viot particle system will be well defined at all time t > 0. The condition τ ∞ = +∞ is clearly fulfilled for continuous time Markov chains with bounded jump rates. In the diffusion process case, criteria have been provided in [8], [33], [65] and [66].
In that case, denote by µ N t the empirical distribution of (Z 1 , · · · ,Z N ) at time t: The following result is obtained in [66] by martingale method.
Theorem 34. Assume that for all N ≥ 2, (Z 1 , · · · , Z N ) is well defined at any time t ≥ 0. Then, for any time t > 0, the sequence of empirical distributions (µ N t ) converges in law to the conditioned distribution P Z 0 (Z t ∈ ·|t < T 0 ), when N goes to infinity.
If moreover (Z 1 , · · · , Z N ) is ergodic, we denote by M N its stationary distribution and by X N its empirical stationary distribution, which is defined by (z 1 , · · · ,z N ) ∈ E * is a random vector distributed with respect to M N . In particular, µ N t converges in law to X N when t → ∞. We refer to [65] for the proof of the following theorem.
Assume moreover that (Z 1 , · · · ,Z N ) is ergodic and that the family of laws of (X N ) N ≥2 is uniformly tight. Then the sequence of random probability measures (X N ) converges weakly to α.
If E is a bounded subset of R d , d ≥ 1, and if Z is a drifted Brownian motion with bounded drift which is killed at the boundaries of E, then the assumptions of Theorems 34 and 35 are fulfilled (see [65]). The proofs of Theorems 34 and 35 are based on a coupling argument. More general (but longer) proofs can also be found in [33] or [66]. In particular, these results provide us a numerical approximation method of the Yaglom limit for such processes.
Let us now consider the Kolmogorov diffusion process X defined in (39). In that case, the existence of the Fleming-Viot particle system remains an open problem because of the unboundedness of the drift coefficient. In order to avoid this difficulty, we introduce the law P ε of the diffusion process with bounded coefficients defined by dX ε t = dB t − q(X ε t )dt ; X 0 ∈ (ε, 1/ε), killed when it hits ε or 1 ε . One can easily show that at any time t ≥ 0, the conditioned distribution of X ε converges to the one of X: P ε X ε t ∈ ·|t < T ε ∧ T 1/ε − − → ε→0 P (X t ∈ ·|t < T 0 ) .
The existence of the Yaglom limit denoted by α and the uniqueness of the QSD for P ε are obtained from Pinsky [49]. The following approximation result is proved in [65] using a compactness-uniqueness argument.
For all N ≥ 2, we denote by (X ε,1 , · · · , X ε,N ) the interacting particle system built as above, with the law P . Since the diffusion process X is a drifted Brownian motion with bounded drift evolving in the bounded interval ] ,1/ [, the interacting particle system (X ε,1 , · · · , X ε,N ) fulfills the assumptions of Theorems 34 and 35. Denoting by µ ,N the empirical distribution of the simulable particle system (X ε,1 , · · · , X ε,N ), we get Then, choosing small enough and N big enough, we get a numerical approximation method for the conditioned distribution and the Yaglom limit of X.
Example 6. Let us now develop this simulation method in the case of the Wright-Fisher diffusion conditioned to be absorbed at 0, which evolves in [0,1[ and is defined by This is a model for a bi-type population in which the second type cannot disappear. In that model, Z t is the proportion of the first type in the population at time t ≥ 0 and 1 − Z t the proportion of the other one. The existence of a Yaglom limit for this process has been proved by Huillet in [35], which also proved that it has the density 2 − 2x with respect to the Lebesgue measure.
Using the approximation method described above with = 0.001 and N = 10000, we obtain numerically the density of the Yaglom limit for Z represented in Figure 13, which is very close to the function x → 2 − 2x and shows the efficiency of the method.