Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift

We study the existence and the exponential ergodicity of a general interacting particle system, whose components are driven by independent diffusion processes with values in an open subset of $\mathds{R}^d$, $d\geq 1$. The interaction occurs when a particle hits the boundary: it jumps to a position chosen with respect to a probability measure depending on the position of the whole system. Then we study the behavior of such a system when the number of particles goes to infinity. This leads us to an approximation method for the Yaglom limit of multi-dimensional diffusion processes with unbounded drift defined on an unbounded open set. While most of known results on such limits are obtained by spectral theory arguments and are concerned with existence and uniqueness problems, our approximation method allows us to get numerical values of quasi-stationary distributions, which find applications to many disciplines. We end the paper with numerical illustrations of our approximation method for stochastic processes related to biological population models.


Introduction
Let D ⊂ R d be an open set with a regular boundary (see Hypothesis 1). The first part of this paper is devoted to the study of interacting particle systems (X 1 ,...,X N ), whose components X i evolve in D as diffusion processes and jump when they hit the boundary ∂D. More precisely, let N ≥ 2 be the number of particles in our system. Let us consider N independent d-dimensional Brownian motions B 1 ,...,B N and a jump measure J (N ) : ∂(D N ) → M 1 (D N ), where M 1 (D N ) denotes the set of probability measures on D N . We build the interacting particle system (X 1 ,...,X N ) with values in D N as follows. At the beginning, the particles X i evolve as independent diffusion processes with values in D defined by dX (i) where q (N ) i is locally Lipschitz on D, such that the diffusion process doesn't explode in finite time. When a particle hits the boundary, say at time τ 1 , it jumps to a position chosen with respect to J (N ) (X 1 τ 1 -,...,X N τn-). Then the particles evolve independently with respect to (1) until one of them hits the boundary and so on. In the whole study, we require the jumping particle to be attracted away from the boundary by the other ones during the jump (in the sense of Hypothesis 2 on J (N ) in Section 2.2). We emphasize the fact that the diffusion processes which drive the particles between the jumps can depend on the particles and their coefficients aren't necessarily bounded (see Hypothesis 1). This construction is a generalization of the Fleming-Viot type model introduced in [5] for Brownian particles and in [20] for diffusion particles. Diffusions with jumps from the boundary have also been studied in [3], with a continuity condition on J (N ) that isn't required in our case, and in [19], where fine properties of a Brownian motion with rebirth have been established.
In a first step, we show that the interacting particle system is well defined, which means that accumulation of jumps doesn't occur before the interacting particles system goes to infinity. Under additional conditions on q (N ) i and D, we prove that the interacting particle system doesn't reach infinity in finite time almost surely. In a second step, we give suitable conditions ensuring the system to be exponentially ergodic. The whole study is made possible thanks to a coupling between (X 1 ,...,X N ) and a system of N independent 1-dimensional reflected diffusion processes. The coupling is built in Section 2.3.
Assume that D is bounded. For all N ≥ 2, let J (N ) be a jump measure and (q (N ) i ) 1≤i≤N a family of drifts. Assume that the conditions for existence and ergodicity of the interacting process are fulfilled for all N ≥ 2. Let M N be its stationary distribution. We denote by X N the associated empirical stationary distribution, which is defined by where (x 1 ,...,x N ) ∈ D N is distributed following M N . Under some bound assumptions on (q (N ) i ) 1≤i≤N,2≤N (see Hypothesis 4), we prove in Section 2.4 that the family of random measures X N is uniformly tight.
In Section 3, we study a particular case: q (N ) i = q doesn't depend on i,N and It means that at each jump time, the jumping particle is sent to the position of a particle chosen uniformly between the N −1 remaining ones. In this situation, we identify the limit of the family of empirical stationary distributions (X N ) N ≥2 . This leads us to an approximation method of limiting conditional distributions of diffusion processes absorbed at the boundary of an open set of R d , studied by Cattiaux and Méléard in [7] and defined as follows. Let U ∞ ⊂ R d be an open set and P ∞ be the law of the diffusion process defined by the SDE and absorbed at the boundary ∂U ∞ . Here B is a d-dimensional Brownian motion and V ∈ C 2 (U ∞ ,R). We denote by τ ∂ the absorption time of the diffusion process (3). As proved in [7], the limiting conditional distribution exists and doesn't depend on x ∈ U ∞ , under suitable conditions which allow the drift ∇V and the set U ∞ to not fulfill the conditions of Section 2 (see Hypothesis 5 in Section 3). This probability is called the Yaglom limit associated with P ∞ . It is a quasi-stationary distribution for the diffusion process (3), which means that P ∞ ν∞ (X ∞ t ∈ dx|t < τ ∂ ) = ν ∞ for all t ≥ 0. We refer to [6,23,25] and references therein for existence or uniqueness results on quasi-stationary distributions in other settings.
Yaglom limits are an important tool in the theory of Markov processes with absorbing states, which are commonly used in stochastic models of biological populations, epidemics, chemical reactions and market dynamics (see the bibliography [29,Applications]). Indeed, while the long time behavior of a recurrent Markov process is well described by its stationary distribution, the stationary distribution of an absorbed Markov process is concentrated on the absorbing states, which is of poor interest. In contrast, the limiting distribution of the process conditioned to not being absorbed when it is observed can explain some complex behavior, as the mortality plateau at advanced ages (see [1] and [32]), which leads to new applications of Markov processes with absorbing states in biology (see [24]). As stressed in [28], such distributions are in most cases not explicitly computable. In [7], the existence of the Yaglom limit is proved by spectral theory arguments, which doesn't allow us to get its explicit value. The main motivation of Section 3 is to prove an approximation method of ν ∞ , even when the drift ∇V and the domain U ∞ don't fulfill the conditions of Section 2.
The approximation method is based on a sequence of interacting particle systems defined with the jump measures (2), for all N ≥ 2. In the case of a Brownian motion absorbed at the boundary of a bounded open set (i.e. q = 0), Burdzy et al. conjectured in [4] that the unique limiting measure of the sequence (X N ) N ∈N is the Yaglom limit ν ∞ . This has been confirmed in the Brownian motion case (see [5], [18] and [26]) and proved in [16] for some Markov processes defined on discrete spaces. New difficulties arise from our case. For instance, the interacting particle process introduced above isn't necessarily well defined, since it doesn't fulfill the conditions of Section 2. To avoid this difficulty, we introduce a cut-off of U ∞ near its boundary. More precisely, let (U m ) m≥0 be an increasing family of regular bounded subsets of U ∞ , such that ∇V is bounded on each U ∞ and such that U ∞ = m≥0 U ∞ . We define an interacting particle process (X m,1 ,...,X m,N ) on each subset U N m , by setting q (1). For all m ≥ 0 and N ≥ 2, (X m,1 ,...,X m,N ) is well defined and exponentially ergodic. Denoting by X m,N its empirical stationary distribution, we prove that We conclude in Section 3.3 with some numerical illustrations of our method applied to the 1-dimensional Wright-Fisher diffusion conditioned to be absorbed at 0, to the Logistic Feller diffusion and to the 2-dimensional stochastic Lotka-Volterra diffusion. 2 A general interacting particle process with jumps from the boundary

Construction of the interacting process
Let D be an open subset of R d , d ≥ 1. Let N ≥ 2 be fixed. For all i ∈ {1,...,N}, we denote by P i the law of the diffusion process X (i) , which is defined on D by and is absorbed at the boundary ∂D. Here B 1 ,...,B N are N independent d-dimensional Brownian motions and q is locally Lipschitz. We assume that the process is absorbed in finite time almost surely and that it doesn't explode to infinity in finite time almost surely.
The infinitesimal generator associated with the diffusion process (5) will be denoted by L We define a system of particles (X 1 ,...,X N ) with values in D N , which is càdlàg and whose components jump from i D i . Between the jumps, each particle evolves independently of the other ones with respect to P i .
Let J (N ) : N i=0 D i → M 1 (D) be the jump measure, which associates a probability measure J (N ) (x 1 ,...,x N ) on D to each point (x 1 ,...,x N ) ∈ N i=1 D i . Let (X 1 0 ,...,X N 0 ) ∈ D N be the starting point of the interacting particle process (X 1 ,...,X N ), which is built as follows: • Each particle evolves following the SDE (5) independently of the other ones, until one particle, say X i 1 , hits the boundary at a time which is denoted by τ 1 . On the one hand, we have τ 1 > 0 almost surely, because each particle starts in D. On the other hand, the particle which hits the boundary at time τ 1 is unique, because the particles evolves as independent Itô's diffusion processes in D. It follows that (X 1 τ 1 -,...,X N τ 1 -) belongs to D i 1 .
• The position of X i 1 at time τ 1 is then chosen with respect to the probability measure J (N ) (X 1 τ 1 -,...,X N τ 1 -). • At time τ 1 and after proceeding to the jump, all the particles are in D. Then the particles evolve with respect to (5) and independently of each other, until one of them, say X i 2 , hits the boundary, at a time which is denoted by τ 2 . As above, we have τ 1 < τ 2 and (X 1 • The position of X i 2 at time τ 2 is then chosen with respect to the probability measure J (N ) (X 1 τ 2 -,...,X N τ 2 -). • Then the particles evolve with law P i and independently of each other, and so on.
The law of the interacting particle process with initial distribution m ∈ M 1 (D N ) will be denoted by P N m , or by P N x if m = δ x , with x ∈ D N . The associated expectation will be denoted by E N m , or by E x if m = δ x . For all β > 0, we denote by S β = inf{t ≥ 0, (X 1 ,...,X N ) ≥ β} the first exit time from {x ∈ D N , x < β}. We set S ∞ = lim β→∞ S β .
The sequence of successive jumping particles is denoted by (i n ) n≥1 , and denotes the strictly increasing sequence of jumping times (which is well defined for all n ≥ 0 since the process is supposed to be absorbed in finite time almost surely). Thanks to the non-explosion assumption on each P i , we have τ n < S ∞ for all n ≥ 1 almost surely. We set τ ∞ = lim n→∞ τ n ≤ S ∞ . The process described above isn't necessarily well defined for all t ∈ [0,S ∞ [, and we need more assumptions on D and on the jump measure J (N ) to conclude that τ ∞ = S ∞ almost surely.
In the sequel, we denote by φ D the Euclidean distance to the boundary ∂D: For all r > 0, we define the collection of open subsets D r = {x ∈ D, φ D (x) > r}. For all β > 0, we set B β = {x ∈ D, x < β}.

Hypothesis 1.
There exists a neighborhood U of ∂D such that 1. the distance φ D is of class C 2 on U, 2. for all β > 0, inf In particular, Hypothesis 1 implies Remark 1. For example, the first part of Hypothesis 1 is fulfilled if D is an open set whose boundary is of class C 2 (see [12,Theorem 4.3]). It is also satisfied by the rectangle with rounded corner defined in Section 3.3.3.
The following assumption ensures that the jumping particle is attracted away from the boundary by the other ones.

Hypothesis 2.
There exists a non-decreasing continuous function f (N ) : R + → R + vanishing at 0 and strictly increasing in a neighborhood of 0 such that, ∀i ∈ {1,...,N}, is a kind of distance from the boundary and we assume that at each jump time τ n , the probability of the event "the jump position X in τn is chosen farther from the boundary than at least one another particle" is bounded below by a positive constant p In that case, the particle on the boundary jumps to one of the other ones, with positive weights. It yields that Hypothesis 2 is fulfilled with p In Section 3, we will focus on the particular case That will lead us to an approximation method of the Yaglom limit (4). Finally, given a jump measure J (N ) satisfying Hypothesis 2 (with p Finally, we give a condition which ensures the exponential ergodicity of the process. In particular, this condition is satisfied if D is bounded and fulfills Hypothesis 1.
If Hypotheses 2 and 3 are fulfilled, then the process (X 1 ,...,X N ) is exponentially ergodic, which means that there exists a probability measure M N on D N such that, where C (N ) (x) is finite, ρ (N ) < 1 and ||.|| T V is the total variation norm. In particular, M N is a stationary measure for the process (X 1 ,...,X N ).
The main tool of the proof is a coupling between (X 1 t ,..., and each i ∈ {1,...,N}. We build this coupling in Subsection 2.2 and we conclude the proof of Theorem 2.1 in Subsection 2.3 . In Subsection 2.4, we assume that D is bounded and that, for all N ≥ 2, we're given J (N ) and a family of drifts (q (N ) i ) 1≤i≤N , such that Hypotheses 1, 2 and 3 are fulfilled. Moreover, we assume that α in Hypothesis 3 doesn't depend on N. Under some suitable bounds on the family (q (N ) i ) 1≤i≤N, N ≥2 , we prove that the family of empirical distributions (X N ) N ≥2 is uniformly tight. It means that, ∀ǫ ≥ 0, there exists a compact set K ⊂ D such that E(X N (D \ K)) ≤ ǫ for all N ≥ 2. In particular, this implies that (X N ) N ≥2 is weakly compact, thanks to [22]. Let us recall that a sequence of random measures (γ N ) N on D converges weakly to a random measure γ on D, if E(γ N (f )) converges to E(γ(f )) for all continuous bounded functions f : D → R. This property will be crucial in Section 3.
We define a sequence of stopping times (θ i n ) n such that More precisely, we set (see Figure  2) and, for n ≥ 1, The sequence (θ i n ) is non-decreasing and goes to τ ∞ ∧ S β almost surely. Let γ i be a 1-dimensional Brownian motion independent of the process (X 1 ,...,X N ) and of the Brownian motion (B 1 ,...,B N ). We set and, for all n ≥ 0, ∇φ D (X i s-) · dB i s has the law of a Brownian motion between times θ i 2n and θ i 2n+1 , thanks to (6). The process ..,N}. Thanks to Hypothesis 2, there exists Q Let us prove that the reflected diffusion process Y β,i defined by (7) fulfills inequality (8) for and we work conditionally to ζ < τ ∞ ∧ S β . By right continuity of the two processes, One can find a stopping time ζ ′ ∈]ζ,τ ∞ ∧ S β [, such that X i doesn't jump between ζ and ζ ′ and such that Y β,i Thanks to the regularity of φ D on B β \ D 2a , we can apply Itô's formula to (φ D (X i t )) t∈[ζ,ζ ′ ] , and we get, for all stopping time t ∈ [ζ,ζ ′ ], But ζ and ζ ′ lie between an entry time of X i to B β \ D a and the following entry time to D 2a .
where the second inequality comes from the positivity of the jumps of φ D (X i ) and from the left continuity of Y β,i , while the third inequality is due to the definition of ζ. Then φ D (X i ) − Y β,i stays non-negative between times ζ and ζ ′ , what contradicts the definition of ζ. Finally, ζ = τ ∞ ∧ S β almost surely, which means that the coupling inequality (8) remains true for all t ∈ [0,τ ∞ ∧ S β [.

Proof of Theorem 2.1
Proof that (X 1 ,...,X N ) is well defined under Hypotheses 1 and 2. Let N ≥ 2 be the size of the interacting particle system and fix arbitrarily its starting point x ∈ D N . Thanks to the non explosiveness of each diffusion process P i , the interacting particle process can't escape to infinity in finite time after a finite number of jumps. It yields that τ ∞ ≤ S ∞ almost surely. Fix β > 0 such that x ∈ B β and define the event C β = {τ ∞ < S β }. Assume that C β occurs with positive probability. Conditionally to C β , the total number of jumps is equal to +∞ before the finite time τ ∞ . There is a finite number of particles, then at least one particle makes an infinite number of jumps before τ ∞ . We denote it by i 0 (which is a random index).
For each jumping time τ n , we denote by σ i 0 n the next jumping time of ) is a continuous diffusion process with bounded coefficients between τ n and σ i 0 n -, then Since the process φ D (X i 0 ) is continuous between τ n and σ i 0 n −, we conclude that φ D (X i 0 τn ) doesn't lie above the support of f , for n big enough almost surely. But the support of f can be chosen arbitrarily close to 0, it yields that φ D (X i 0 τn ) goes to 0 almost surely conditionally to C β .
Let us denote by (τ i 0 n ) n the sequence of jumping times of the particle i 0 . We denote by A n the event where f (N ) is the function of Hypothesis 2 . We have, for all 1 ≤ k ≤ l, where, by definition of the jump mechanism of the interacting particle system, by Hypothesis 2. By induction on l, we get It means that, for infinitely many jumps τ n almost surely, one can find a particle j such that Because there is only a finite number of other particles, one can find a particle, say j 0 (which is a random variable), such that , for infinitely many n ≥ 1.
In particular, lim n→∞ φ D (X i 0 τn ),f (N ) (φ D (X j 0 τn )) = (0,0) almost surely. But (f (N ) ) −1 is well defined and continuous near 0, then Using the coupling inequality of Proposition 2.2, we deduce that Then, conditionally to C β , Y β,i 0 and Y β,j 0 are independent reflected diffusion processes with bounded drift, which hit 0 at the same time. This occurs for two independent reflected Brownian motions with probability 0, and then for Y β,i 0 and Y β,j 0 too, by the Girsanov's Theorem. That implies P x (C β ) = 0.
If the first part of Hypothesis 3 is fulfilled, one can defined the coupled reflected diffusion Y ∞,i , which fulfills inequality (8) with a = α and for all t ∈ [0,τ ∞ ∧ S ∞ [= [0,τ ∞ [. Then the same proof leads to Finally, we deduce that τ ∞ = ∞ almost surely.
where φ is a regular function. In our case of a drifted Brownian motion, φ is equal to 1 and Y i is a reflected drifted Brownian motion independent of the others particles. But in the general case, the Y i are general orthogonal semi-martingales. It yields that the generalization of the previous proof reduces to the following hard problem (see [31, Question 2, page 217] and references therein): "Which are the two-dimensional continuous semi-martingales for which the one point sets are polar ?". Since this question has no general answer, it seems that the previous proof doesn't generalize immediately to general uniformly elliptic diffusion processes.
We emphasize the fact that the proof of the exponential ergodicity can be generalized (as soon as τ ∞ = S ∞ = +∞ is proved), using the fact that (Y 1 t ,...,Y N t ) t≥0 is a time changed Brownian motion with drift and reflection (see [31, Theorem 1.9 (Knight)]). This time change argument has been developed in [20], with a different coupling construction. This change of time can also be used in order to generalize Theorem 2.3 below, as soon as the exponential ergodicity is proved.
Proof of the exponential ergodicity. It is sufficient to prove that there exists n ≥ 1, ǫ > 0 and a non-trivial probability ϑ on D N such that are defined in Hypothesis 3, and such that where κ is a positive constant and τ ′ = min{n ≥ 1, (X 1 ) n∈N is aperiodic (which is obvious in our case) and fulfills (9) and (10), then it is geometrically ergodic. But, thanks to [13, Theorem 5.3 p.1681], the geometric ergodicity of this Markov chain is a sufficient condition for (X 1 ,...,X N ) to be exponentially ergodic.
We assume without loss of generality that K (N ) 0 ⊂ D α/2 (where α is defined in Hypothesis 3). Let us set . Thanks to Hypothesis 3, ϑ is a non-trivial probability measure. Moreover, (9) is clearly fulfilled with n = 1 and ǫ = N i=1 inf x∈Dα P i (X > 0 is defined in Hypothesis 3. It yields that It yields that there exists κ > 0 such that (10) is fulfilled.

Uniform tightness of the empirical stationary distributions
In this part, the open set D is supposed to be bounded. Assume that a jump measure J (N ) and a family of drifts (q (N ) i ) i=1,...,N are given for each N ≥ 2.
Hypothesis 4. Hypotheses 1 and 2 are fulfilled for each N ≥ 2 and Hypothesis 3 is fulfilled with the same α for each N ≥ 2. Moreover, there exists r > 1 such that For all N ≥ 2, we denote by m N ∈ M 1 (D N ) the initial distribution and by µ N (t,dx) the empirical distribution of the N-particles process defined by the jump measure J (N ) and the family (q (N ) i ) i∈{1,...,N } . Its stationary distribution is denoted by M N and its empirical stationary distribution is denoted by X N : where (x 1 ,...,x N ) is a random vector in D N distributed following M N .
Thanks to the Girsanov's Theorem, we have where (w 1 ,...,w N ) is a N-dimensional Brownian motion. By the Cauchy Schwartz inequality, we get where the second inequality occurs, since 0 ≤ δ w i i w i t , whose expectation is 1. Taking the expectation in (11), it yields that Thanks to Hypothesis 4, there exists goes to 0 when r → 0, so that the family of random measures (µ N (t,dx)) N ≥2 is uniformly tight.
If we set m N equal to the stationary distribution M N , then we get by stationarity that X N is distributed as µ N (t,.), for all N ≥ 2 and t > 0. Finally, the family of empirical stationary distributions (X N ) N ≥2 is uniformly tight.

Yaglom limit's approximation
We consider now the particular case J (N ) (x 1 ,...,x N ) = 1 N −1 N k=1,k =i δ x k : at each jump time, the particle which hits the boundary jumps to the position of a particle chosen uniformly between the N − 1 remaining ones. We assume moreover that q (N ) i = q doesn't depend on i,N. In this framework, we are able to identify the limiting distribution of the empirical stationary distribution sequence, when the number of particles tends to infinity. This leads us to an approximation method of the Yaglom limits (4), including cases where the drift of the diffusion process isn't bounded and where the boundary is neither regular nor bounded.
Let U ∞ be an open domain of R d , with d ≥ 1. We denote by P ∞ the law of the diffusion process defined on U ∞ by 12) and absorbed at the boundary ∂U ∞ . Here B is a d-dimensional Brownian motion and V ∈ C 2 (U ∞ ,R). We assume that Hypothesis 5 below is fulfilled, so that the Yaglom limit exists and doesn't depend on x, as proved by Cattiaux and Méléard in [7,Theorem B.2]. We emphasize the fact that this hypothesis allows the drift ∇V of the diffusion process (12) to be unbounded and the boundary ∂U ∞ to be neither of class C 2 nor bounded. In particular, the results of the previous section aren't available in all generality for diffusion processes with law P ∞ .
Hypothesis 5. We assume that 5. There exists R 0 > 0 such that Here p U∞ 1 is the transition density of the diffusion process (12) with respect to the Lebesgue measure.
According to [7], the second point implies that the semi-group induced by P ∞ is ultracontractive. The assumptions 1-4 imply that the generator associated with P ∞ has a purely discrete spectrum and that its minimal eigenvalue −λ ∞ is simple and negative. The last assumption ensures that the eigenfunction associated with −λ ∞ is integrable with respect to e −2V (x) dx. Finally, Hypothesis 5 is sufficient for the existence of the Yaglom limit (13).

Remark 5.
For example, it is proved in [7] that Hypothesis 5 is fulfilled by the Lotka-Volterra system studied numerically in Subsection 3.3.3. Up to a change of variable, this system is defined by the diffusion process with values in U ∞ = R 2 + , which satisfies (14) and is absorbed at ∂U ∞ . Here B 1 ,B 2 are two independent one-dimensional Brownian motions and the parameters of the diffusion process fulfill condition (30).
In order to define the interacting particle process of the previous section, we work with diffusion processes defined on U m , m ≥ 0. More precisely, for all m ≥ 0, we denote by P m the law of the diffusion process defined on U m by 15) and absorbed at the boundary ∂U m . Here B is a d-dimensional Brownian motion and q m : U m → R is a continuous function. We denote by L m the infinitesimal generator of the diffusion process with law P m . For all m ≥ 0, the diffusion process with law P m clearly fulfills the conditions of Section 2. For all N ≥ 2, let (X m,1 ,...,X m,N ) be the interacting particle process defined by the law P m between the jumps and by the jump measure By Theorem 2.1, this process is well defined and exponentially ergodic.
For all m ≥ 0 and all N ≥ 2, we denote by µ m,N (t,dx) the empirical distribution of (X m,1 t ,...,X m,N t ), by M m,N the stationary distribution of (X m,1 ,...,X m,N ) and by X m,N the associated empirical stationary distribution.
We are now able to state the main result of this section.
in the weak topology of random measures, which means that, for all bounded continuous function f : In Section 3.1, we fix m ≥ 0 and we prove that the sequence (X m,N ) N ≥2 converges to a deterministic probability ν m when N goes to infinity. In particular, we prove that ν m is the Yaglom limit associated with P m , which exists by [7]. In Section 3.2, we conclude the proof, proceeding by a compactness/uniqueness argument: we prove that (ν m ) m≥0 is a uniformly tight family and we show that each limiting probability of the family (ν m ) m≥0 is equal to the Yaglom limit ν ∞ . The last Section 3.3 is devoted to numerical illustrations of Theorem 3.1.
3.1 Convergence of (X m,N ) N ≥2 , when m ≥ 0 is fixed Proposition 3.2. Let m ≥ 0 be fixed and let q m : U m → R be a continuous function. Assume that µ m,N (0,dx) converges in the weak topology of random measure to a random probability measure µ m with values in M 1 (U m ), when N → ∞. Then, for all T ≥ 0, µ m,N (T,dx) converges in the weak topology of random measure to P m µm (X T ∈ .|X T ∈ U m ) when N goes to infinity.
Moreover, if there exists ν m ∈ M 1 (U m ) such that then the sequence of empirical stationary distributions (X m,N ) N ≥2 converges to ν m in the weak topology of random measures when N goes to infinity.
Remark 6. In Proposition 3.2, ν m is the Yaglom limit and the unique quasi-stationary distribution associated with P m . For instance, each of the two following conditions is sufficient for the existence of such a measure: 1. If q m = 1 Um ∇V , by [7]. This is the case of Theorem 3.1.
Proof of Proposition 3.2. We set where A N t = ∞ n=1 1 τn≤t denotes the number of jumps before time t. Intuitively, we introduce a loss of 1/N of the total mass at each jump, in order to approximate the distribution of the diffusion process (15) without conditioning. We will come back to the study of µ m,N and the conditioned diffusion process by normalizing ν m,N .
From the Itô's formula applied to the semi-martingale µ m,N (t,ψ) = 1 where M c,N (t,ψ) is the continuous martingale Applying the Itô's formula to the semi-martingale ν m,N (t,ψ), we deduce from (17) that Where we have That implies It yields that, for all smooth functions Ψ(t,x) vanishing at the boundary of U m , Let T > 0 be fixed. For all δ > 0, define Ψ δ (t,x) = P m T −t P m δ f (x), where f ∈ C 2 (U m ) and (P m s ) s≥0 is the semigroup associated with P m : P m s f (x) = E x (f (X Um s )). Then Ψ δ vanishes on the boundary, is smooth, and fulfills thanks to Kolmogorov's equation (see [14,Proposition 1.5 p.9]). It yields that where c m > 0 is a positive constant. The last inequality comes from [30,Theorem 4.5] on gradient estimates in regular domains of R d . The jumps of the martingale M j,N (t,Ψ δ ) are smaller than 2 N Ψ δ ∞ , then Taking t = T and δ = 1 N , we get from (18), (19) and (20) that Assume that f vanishes at ∂U m , so that f belongs to the domain of L m . Then P m By assumption, the family of random probabilities (ν m,N (0,.)) N ≥2 = (µ m,N (0,.)) N ≥2 converges to µ m . We deduce from (21) that for all f ∈ C 2 (U m ) vanishing at boundary. But the family ν m,N (T,.) N ≥2 is uniformly tight by Theorem 2.3 . It yields from (22) that its unique limiting distribution is µ m (P m T .). In particular, ) . But µ m (P m T 1 Um ) never vanishes almost surely, so that The family of random probabilities (X m,N ) N ≥0 is uniformly tight, by Theorem 2.3. Let X m be one of its limiting probabilities. By definition, there exists a strictly increasing map ϕ : N → N, such that X m,ϕ(N ) converges in distribution to X m when N → ∞. By stationarity, X m,ϕ(N ) has the same law as µ m,ϕ(N ) (T,.), which converges in distribution to converges almost surely to ν m when T → ∞, by (16). We deduce from this that X m has the same law as ν m . As a consequence, the unique limiting probability of the uniformly tight family (X m ) N is ν m , which allows us to conclude the proof of Proposition 3.2.

Convergence of the family (ν m ) m≥0
Proposition 3.3. Assume that Hypothesis 5 is fulfilled and that q m = ∇V 1 Um . Then the sequence (ν m ) m≥0 converges weakly to the Yaglom limit ν ∞ when m → ∞.
Remark 7. Since q m = ∇V 1 Um , the operator L m is symmetric with respect to the measure e −2V (x) dx, but this isn't directly used in the proof of Proposition 3.3. We mainly use inequalities from [7] that are implied by the ultra-contractivity of P ∞ and the third point of Hypothesis 5. However, it seems hard to generalize this last hypothesis and its implications to diffusions with non-gradient drifts.
Proof of Proposition 3.3. For all m ≥ 0 and m = ∞, it has been proved in [7] that −L m * has a simple eigenvalue λ m > 0 with minimal real part, where L m * is the adjoint operator of L m . The corresponding normalized eigenfunction η m is strictly positive on U m , belongs to C 2 (U m ,R) and fulfills where dσ(x) = e −2V (x) dx.
The Yaglom limit ν m is given by In order to prove that (ν m ) m≥0 converges to ν ∞ , we show that (λ m ) m≥0 converges to λ ∞ . Then we prove that (η m 1 Um dσ) m≥0 is uniformly tight. We conclude by proving that every limiting point ηdσ is a nonzero measure proportional to η ∞ dσ. For all m ≥ 0 or m = ∞, the eigenvalue λ m of −L m * is given by (see for instance [34,chapter XI,part 8]) where C ∞ 0 (U m ) is the vector space of infinitely differentiable functions with compact support in U m and f,g σ,m = Um f (u)g(u)dσ(u). For all φ ∈ C ∞ 0 (U ∞ ), the support of φ belongs to U m for m big enough, then Let us show that the family (η m 1 Um dσ) m≥0 is uniformly tight. Fix an arbitrary positive constant ǫ > 0 and let us prove that one can find a compact set K ǫ ⊂ U ∞ which fulfills U∞\Kǫ ǫ m 1 Um dσ ≤ ǫ, ∀m ≥ 0.
which is smaller than ǫ/2 for a good choice of K, say K ′ ǫ , since the integral at the right-hand side is finite by Hypothesis 5. On the other hand where κ = sup m≥0 η m e −V ∞ < ∞ thanks to [7], and λ m ≤ λ ∞ for all m ≥ 0. But the integral on the right-hand side is well defined by Hypothesis 5, then one can find a compact set K ′′ ǫ such that (26) is bounded by ǫ/2. We set K ǫ = K ′ ǫ ∪ K ′′ ǫ so that (25) is fulfilled. Since inequality (25) occurs for all ǫ > 0, the family (η m dσ) m≥0 is uniformly tight. Moreover, η m dσ has a density with respect to the Lebesgue measure, which is bounded by κe −V , uniformly in m ≥ 0. Then it is uniformly bounded on every compact set, so that every limiting distribution is absolutely continuous with respect to the Lebesgue measure.
Let ηdσ be a limiting measure of (η m dσ) m≥0 . For all φ ∈ C ∞ 0 (U ∞ ,R), the support of φ belongs to U m for m big enough, then Thanks to the elliptic regularity Theorem, η is of class C 2 and fulfills L ∞ * η = −λ ∞ η. But the eigenvalue λ ∞ is simple, then η is proportional to η ∞ . Let β ≥ 0 be the non-negative constant such that η = βη ∞ . In particular, there exists an increasing function φ : N → N such that η φ(m) dσ converges weakly to βη ∞ dσ.
Let us prove that β is positive. For all compact subset K ⊂ U ∞ , we have where κ = sup m≥0 η m e −V ∞ < ∞. For all m ≥ 0 and all R > 0, where G and G are defined in Hypothesis 5. Let us prove that η m ,Gη m σ,m is uniformly bounded in m ≥ 0. For all x ∈ U m , (24) leads to Then where the second equality is a consequence of the Green's formula (see [2,Corollary 3.2.4]). But G(R) goes to +∞ when R → ∞, then there exists R 1 > 0 such that 1 G(R 1 ) η m ,1 |x|≥R 1 Gη m σ,m ≤ 1 4 . Since κ = sup m≥0 η m e −V ∞ < ∞, we deduce from (28) that But one can find a compact subset K 1 ⊂ U ∞ such that U∞ 1 {|x|<R 1 }\K 1 dx ≤ 1 4κ 2 , then we have from (27) It yields that β > 0 and Proposition 3.3 follows.

The Wright-Fisher case
The Wright-Fisher with values in ]0,1[ conditioned to be absorbed at 0 is the diffusion process driven by the SDE and absorbed when it hits 0 (1 is never reached). Huillet proved in [21] that the Yaglom limit of this process exists and has the density 2 − 2x with respect to the Lebesgue measure. In order to apply Theorem 3.1, we define P ∞ as the law of X ∞ . = arccos(1 − 2Z . ). Then P ∞ is the law of the diffusion process with values in U ∞ =]0,π[, driven by the SDE dX ∞ t = dB t − 1 − 2 cos X ∞ t 2 sin X ∞ t dt, X ∞ 0 = x ∈]0,π[, absorbed when it hits 0 (π is never reached). One can easily check that this diffusion process fulfills Hypothesis 5. We denote by ν ∞ its Yaglom limit. For all m ≥ 1, we define U m =] 1 m ,π− 1 m [. Let P m and ν m be as in Section 3. We proceed to the numerical simulation of the N-interacting particle system (X m,1 ,...,X m,N ) with m = 1000 and N = 1000. This leads us to the computation of E(X m,N ), which is an approximation of ν ∞ . After the change of variable Z. = 2 cos(X.), we see on Figure 3 that the simulation is very close to the expected result (2 − 2x)dx, which shows the efficiency of the method.

The logistic case
The logistic Feller diffusion with values in ]0, + ∞[ is defined by the stochastic differential equation dZ t = Z t dB t + (rZ t − cZ 2 t )dt, Z 0 = z > 0, and absorbed when it hits 0. Here B is a 1-dimensional Brownian motion and r,c are two positive constants. In order to use Theorem 3.1, we make the change of variable X. = 2 √ Z..
We denote by P ∞ its law. Cattiaux et al. proved in [6] that Hypothesis 5 is fulfilled in this case. Then the Yaglom limit ν ∞ associated with P ∞ exists and one can apply Theorem 3.1 with U m =] 1 m ,m[ for all m ≥ 1. For all N ≥ 2, we denote by P m the law of the diffusion process restricted to U m and by X m,N the empirical stationary distribution of the N-interacting particle process associated with P m .
We've proceeded to the numerical simulation of the interacting particle process for a large number of particles and a large value of m. This allows us to compute E(X m,N ), which gives us a numerical approximation of ν ∞ , thanks to Theorem 3.1.
In the numerical simulations below, we set m = 10000 and N = 10000. We compute E(X m,N ) for different values of the parameters r and c in (29). The results are graphically represented in Figure 4. As it could be wanted for, greater is c, closer is the support of the QSD to 0. We thus numerically describe the impact of the linear and quadratic terms on the Yaglom limit.

Stochastic Lotka-Volterra Model
We apply our results to the stochastic Lotka-Volterra system with values in D = R 2 + studied in [7], which is defined by the following stochastic differential system where (B 1 ,B 2 ) is a bi-dimensional Brownian motion. We are interested in the process absorbed at ∂D.
In [7], this case was called the weak cooperative case and the authors proved that it is a sufficient condition for Hypothesis 5 to be fulfilled. Then the Yaglom limit ν ∞ = lim t→+∞ P ∞ x (X ∞ ∈ .|t < τ ∂ ) is well defined and we are allowed to apply Theorem 3.1. For each m ≥ 1, we define U m as it is described on Figure 5. With this definition, it is clear that all conditions of Theorems 2.1 and 3.1 are fulfilled.
We choose m = 10000 and we simulate the long time behavior of the interacting particle process with N = 10000 particles for different values of c 12 and c 21 . The values of the other parameters are r 1 = 1 = r 2 = 1, c 11 = c 22 = 1, γ 1 = γ 2 = 1. The results are illustrated on Figure 6. One can observe that a greater value of the cooperating coefficients −c 12 = −c 21 leads to a Yaglom limit whose support is further from the boundary and covers a smaller area. In other words, the more the two populations cooperate, the bigger the surviving populations are.