Recent advances in the long-time analysis of killed degenerate processes and their particle approximation

We review some recent results of quantitative long-time convergence for the law of a killed Markov process conditioned to survival toward a quasi-stationary distribution, and on the analogous question for the particle systems used in practice to sample these distributions. With respect to the existing literature, one of the novelties of these works is the degeneracy of the underlying process with respect to classical elliptic diffusion, namely it can be a non-elliptic hypoelliptic diffusion, a piecewise deterministic Markov process or an Euler numerical scheme.


Killed Markov processes
A killed process on a space D is a Markov process (X t ) t⩾0 on D ∪{∂} where ∂ / ∈ D is an arbitrary cemetery point which is an absorbing point, i.e.X t = ∂ for all t ⩾ τ = inf{s ⩾ 0, X s = ∂}.In many cases of interest, it is obtained from a Markov process (Y t ) t⩾0 on a state S by setting X t = Y t for t < τ and X t = ∂ for t ⩾ τ where τ is defined either by a so-called soft-killing mechanism, namely, given E a standard exponential random variable independent from Y , τ = inf t ⩾ 0, E ⩽ t 0 λ(Y s )ds , with a death rate λ : S → R + , in which case D = S, or a so-called hard-killing mechanism, namely where D is a given domain of S. Since we are only interested in the process before absorption, in these cases it is equivalent to work with either X or Y .Some objects of interest are then, for all t ⩾ 0, p t = P(τ ≤ t) the extinction probability, η t the law of the killed process and µ t the law of the process conditioned to survival, given by for all measurable sets A of D.
In most cases, extinction is almost sure, i.e. τ < +∞ almost surely.In that case, the only invariant measure of X is the Dirac mass at ∂, which is not interesting.In fact, for killed processes, the relevant analogous to invariant measures is given by the quasi-stationary distributions (QSD), which by definitions are the probability measures ν on D such that Similarly, the classical question of convergence toward equilibrium for Markov processes is replaced, for killed process, by the question of the limit of µ t toward a QSD as t → +∞ when the initial distribution is not a QSD.
An important property is that, starting from a QSD, the extinction time follows an exponential distribution.In other words, if ν is a QSD, then there exists θ > 0 such that, if X 0 ∼ ν, then p t = 1 − e −θt for all t ⩾ 0. Equivalently, η t = e −θt ν for all t ⩾ 0, which means that ν is an eigenvector of the killed semi-group.In fact, when D is a finite set, the existence of a QSD stems from the Perron-Frobenius theorem.
For general and extensive references on QSD, we refer the interested reader to the book [22] or the survey [55].

Some motivations
Among other things, the study of killed processes and quasi-stationary distributions is related to metastability phenomena, i.e. cases where a Markov process stays for very long times in some transient states before reaching stationarity.The convergence to a QSD within a transient state is fast and observed before the slow transitions.This arises in particular in population models and stochastic algorithms in molecular dynamics, as we describe in the following.
In population dynamics, a Markov process (X t ) t≥0 is generally used to model the number of individuals.This can for example be a diffusion process on a bounded set as in the case of the evolution of gene populations in the continuous Wright-Fischer model [10,Page 368]  by the chemostat model [17,21] where X t = (N t , S t ) ∈ N × R + is a couple of two components: a discrete one N t and a continuous one S t .The discrete component represents the number of bacteria in a vessel and evolves as a pure jump process with the addition or removal of an individual at exponential clocks.Parameters of the jump times depend on the continuous component which evolves according an ordinary differential equation whose parameters depend on the discrete component.See [17,21] for details and Figure 1 for an example of trajectory.This process is not irreducible, belongs to a non-compact space and is hybrid, in the sense that its transition kernel does not admit a density with respect to a simple reference measure.Whatever the parameters, the population of bacteria always hits 0 in finite time (as a consequence of a Borel-Cantelli type arguments since resources are finite), however, when the vessel is large the dynamics is close to a deterministic model having a non-trivial equilibrium.These two antagonistic behaviors, which are common for population models, should be explained by a quasi-stationary behavior.
As we shall see here, the study of QSD can also be useful in molecular dynamics methods for the sampling of the time-evolution of complex particle systems.This sampling is done using numerical schemes which are used in a wide range of application fields (biology, chemistry, materials science for instance).They allow us to compute dynamical quantities of interest for the system without resorting to the actual experiment.For instance, one would like to sample the trajectories of a system which goes from one state to another one, in order to compute the typical time to observe such transitions.This question arises for instance when we want to sample the phase transitions of a thermostated (fixed temperature) molecular system, whose evolution can be modeled by the Langevin dynamics: where (x t , v t ) ∈ R d × R d denotes the positions and momenta of the particles at time t ≥ 0, M is the mass matrix, F is the interaction force (often F is conservative, i.e.F = −∇V for some potential function V ), γ > 0 is a friction parameter, β = (k B T ) −1 is proportional to the inverse temperature and (B t ) t≥0 is a standard Brownian motion.When the friction coefficient γ goes to infinity, it is known that its rescaled position coordinate (x γt ) t≥0 converges to the solution (x t ) t≥0 of the overdamped Langevin Dynamics (see for instance [48, Section 2.2.4]): Phase transitions are known to be often metastable, which means that the timestep used for the sampling of the dynamics is much smaller than the timescale associated with the transition event.As a practical consequence, the number of iterations needed to sample the transition event will be far too large to be done through naive simulation.The metastability comes from the existence of energetic barriers, namely the trajectory to leave the state requires to climb above a saddle point of the potential energy V .
Several methods have been developed in order to circumvent this issue, among which the Parallel Replica method [62].The idea behind Parallel Replica is that, as already mentioned, when the system remains trapped inside a metastable state for a sufficiently long time, its distribution becomes close to the QSD of the state.Furthermore, starting from the QSD, the first exit from the state satisfies some well-known properties which allow for the sampling in parallel of the first exit event, thus reducing the wall clock time simulation in Parallel Replica compared to standard simulation.
The main ingredients for the mathematical justification of Parallel Replica are the existence of a QSD in the state considered, as well as a long time convergence of the dynamics trapped inside the state towards the QSD.Such proofs were carried out in the literature in the case of the overdamped Langevin (4) dynamics in [45] but had yet to be extended to the Langevin dynamics (3).

Past results and limitations
As far as continuous-time processes in a continuous space are concerned, there is a wide literature on quasi-stationary distributions for elliptic diffusions like the overdamped Langevin dynamics (4) on smooth bounded domains.We refer for example to [8,9,11,15,35,45,57]. Let us state a theorem to summarize the main results known for the overdamped Langevin.
The infinitesimal generator of the overdamped Langevin dynamics is given by: with formal adjoint L * in L 2 (dxdv) given by: The following is based on [11,35,41,45].
Theorem 1 (QSD of the overdamped Langevin process).Let Then there exists a unique QSD µ on O for the killed process in O, (X t ) t≥0 , solving (4).Furthermore, (1) there exists ψ ∈ C 2 (O) ∩ C b (O) such that µ(dq) = ψ(q)dq, where dq is the Lebesgue measure on R d , (2) Span(ψ) is the eigenspace associated with the smallest eigenvalue θ of the operator −L * with homogeneous Dirichlet boundary conditions on ∂O, (3) there exist C > 0 and α > 0 such that for all probability measures ν on O, for all t ≥ 0, where µ t is the law of the process (X t ) t≥0 conditioned to survival, as defined in (2), satisfying X 0 ∼ ν and ∥ • ∥ T V is the total-variation norm on the space of bounded signed measures on R d .Now, additional difficulties arise when the process is not an elliptic diffusion process, which is the case of the processes mentioned in the previous section: the kinetic Langevin dynamics ( 3) is an hypoelliptic non-elliptic diffusion, and the chemostat model of [17,21] is a piecewise deterministic Markov process, lacking the regularization properties of diffusions.Moreover, the study of hardkilled process with an unbounded domain (which is the case for the Langevin process since the velocity is unbounded) is also difficult.For instance even for the very simple Ornstein-Uhlenbeck process killed at the origin, there is no uniqueness of the QSD [50].
Multiple approaches were used in the literature to obtain Theorem 1.To name just one, one can apply the Krein-Rutman theorem to the semigroup of (4) killed outside of O, which is compact on L 2 (O), to obtain the QSD.However, such compactness property is not clear once the infinitesimal generator is not elliptic or the domain is unbounded as is the case for (3) on the domain D = O × R d .Two approaches are presented in Sections 3 and 4 to establish this property in cases that cover the kinetic Langevin process (3).Section 2 presents an alternative method, based on a Doeblin-type minorization condition rather than compactness, to get a result similar to the Krein-Rutman theorem in non-compact space.As can be seen by comparing the general results stated in Theorems 2 and 6, the two methods share some similarities (a Lyapunov function is required in both cases) but each have their interest: the conditions of Theorem 6 are easier to check in practice, but the process is required to be Feller, which is not the case for Theorem 2, while in the latter a non-trivial minorization condition has to be established.
Finally, notice that, when it comes to sample the QSD (for instance as the first step of the Parallel Replica algorithm), a naive Monte Carlo method, consisting of running N independent copies of the process up to a long time t and taking the empirical distribution of the processes which have not been killed, may suffer from a high variance since the size of the sample shrinks with times (possibly, all the processes may have already died at time t).For this reason, systems of interacting particle are used in practice: when a process dies, another process is chosen at random uniformly and is duplicated, which maintains the size of the sample, at the cost of a small bias.The particle system is then a non-killed Markov process, which converges in long time to an invariant measure which is close, for large N , to the law of N independent realizations of the QSD.Although many results are available in various frameworks for the convergence of the system of particles toward the killed process as N → ∞, either at a fixed time t or at stationary (see [16,23] and references within) obtaining quantitative convergence as t → ∞ uniformly in N is challenging and cannot be obtained directly from the long-time convergence of the killed process.Such a result, in the case of an elliptic diffusion, is presented in Section 5.The effect of the time-discretization of the continuous-time diffusion, which is a very classical topic for the invariant measure of Markov processes but not for QSD, is also addressed, closing the gap between the theoretical results such as Theorem 1 and the algorithms used in practice to sample the QSD.

2.
Harris-type method for non-conservative processes 2.1.Harris theorem type assumption for Krein-Rutman theorem type conclusion To study (η t ) or (µ t ) defined in (2), it is usual to introduce the family of operators (M t ) defined for all t ≥ 0 by where f is a bounded type function and x ∈ D. It is easy to show that (M t ) defines a positive semigroup, but in contrast with Markov semigroup, it is not conservative in the sense that in general, where 1 is the constant function equal to 1.This particularity is a drawback because there is a rich base of tools to study conservative semigroups as for instance Doeblin contraction or Lyapunov methods.More precisely, a Markov semigroup (P t ) satisfies the Doeblin minorization condition if there exists t 0 , ϵ > 0 and a probability measure such that for all x ∈ X δ x P t0 ≥ ϵν.
This induces uniqueness of an invariant distribution and exponential convergence to it in total variation distance.When the state is not compact Assumption (6) generally does not apply on the whole space.Consequently, it is only supposed to hold on compact sets and an additional assumption is required to prove that the process remains in a compact set for sufficiently long times.This tightness condition is often verified through the existence of a so-called Lyapunov function: there exists a positive function V , such that for a certain t > 0, for some constants C, D > 0 with C < 1.This is often verified by the drift condition: where L is the generator of (P t ) and c, d > 0. With this assumption, Doeblin minorization only need to hold on the sublevel set of V to obtain exponential convergence to a unique invariant measure.
Several works aimed at generalizing these efficient techniques for non-conservative semi-groups.Hilbert metric and Birkhoff contraction yield another powerful method for the analysis of semigroups, which has been well developed [7,56,59].More recently, let us cite for instance [24,43,44] in a general and continuous-time setting or works on discrete-time Feynman-Kac semigroups [23,26].Some works associate Doeblin-Harris techniques with Krein-Rutmann theorem and then often need strong Feller properties [31,36,46]; See Section 3 and 4. Unfortunately, these results do not apply or are difficult to apply for piecewise deterministic dynamics as e.g. the chemostat model described in introduction.
In [2,3,18], some theoretical results were developed to generalize Doeblin-Harris methods for non-conservative semi-groups mainly in case of branching dynamics (i.e. in case where M t 1 > 1).These results are related to [12,13,20,21,25].A generalization of Harris Theorem is for instance described in [3, Theorem 2.1] and is given below.
Theorem 2 (Generalized Harris Theorem).(i) If there exist a couple of positive functions (V, ψ), ) for all positive and measurable function f , then, there exists a unique triplet (γ, h, θ) of eigenelements of M with γ(h) = 1, i.e. satisfying for all t ≥ 0 γM t = e −θt γ and M t h = e −θt h.
By differentiating (8), the triplet (γ, h, θ) is a triplet of eigenelements for the infinitesimal generator A of (M t ) t≥0 , that is γA = −θγ and Ah = −θh where the (unbounded) operator A is defined by The assumptions of Theorem 2 hold on compact space when (M t ) admits a (upper and lower) bounded density with respect to a fixed measure ν.It is known since the 60's, by Birkhoff [7] for instance, that this type of regularity hypothesis implies the existence of eigen-elements and the convergence.Unfortunately, assuming bounded density does not include simple examples such as the chemostat model.
Assumption (A4) overcomes this problem and allows to treat such an example.This type of condition seems to appear for the first time in [25] without the Lyapunov condition.Article [25] only states a contraction inequality, but this is the main step of the proof for such results.In [12], Figure 2. A sample path of the process described in Equation 11.
Champagnat and Villemonais highlight these conditions and prove that they are in fact equivalent to the exponential convergence (9).They deduce the existence of eigenelements in the context of absorbed Markov processes (sub-conservative semigroups).In [2], this theorem is extended for general inhomogeneous-time and non-conservative semi-groups.It includes for instance periodic and asymptotically homogeneous ones.
As explained before, requiring a Doeblin type minorization assumption on the whole space is often too strong when the state space is not compact, Assumptions (A1) and (A2) then relax this constraint.This Lyapunov assumption seems to appear for the first time in [20,Theorem 4.2] and in the proof of [21,Theorem 4.1].However, combining this assumption with the assumptions (A3) and (A4) to obtain uniqueness of eigenelements and exponential convergence comes from [13].This significant result has been demonstrated in the framework of sub-conservative semigroups.Result [3, Theorem 2.1] generalizes this result for non-conservative semi-groups, and non-bounded functions ψ, and prove that these assumptions are in fact necessary.Let us also point out the follow-up paper [14] of [3,13] on these aspects.However, Assumption (A4) of Theorem 2 seems simpler to verify that the analogous of [14] and the proof of [3] leads some quantitative estimates.More precisely, we have that γ(V ) < +∞, γ ≫ ν (in the Radon-Nikodym sense), where c 1 , c 2 > 0, q ∈ (0, 1) and the others constants come from Theorem 2.
To conclude this section let us mention that from [18], Condition (A4) is verified when there exist t 0 , ϵ > 0 such that inf x,y∈K This condition will be verified in almost all examples of Subsection 2.2

Some (very) simple examples
In this section, we detail several examples where one of the assumptions does not hold in order to highlight their role.
Note however that all assumptions are naturally connected and that one can postpone the difficulty of an assumption by changing the compact K or the Lyapunov functions for instance.

Aperiodicity and Assumption (A3)
In addition to the irreducibly assumption, Assumption (A3) hides an aperiodicity assumption.Indeed the (unabsorbed) Markov process (X t ) t≥0 on [0, 1) defined, for all t ≥ 0, by verifies all the assumptions except (A3).Even if the trajectories pass into all the points of [0, 1), it does not pass at the same time t 0 uniformly over its stating point.See Figure 2.
A less simple and more interesting example is given by the growth-fragmentation process studied in [6,34] (in a different formalism).This process grows exponentially at a certain rate α and has multiplicative jumps from x to rx, where r ∈ (0, 1) is a fixed deterministic value, at a rate B(x) depending on its state.We can kill the process at some rate λ as defined in Equation 1.Its generator reads for smooth functions f and x > 0. See Figure 3 In contrast with the previous trivial deterministic process, this growth-fragmentation has a random dynamics but its law is concentrated in fixed time in a discrete space (when the starting distribution is a Dirac mass).It then can not converge to its unique QSD, which has a Lebesgue density.

Probability of survival and Assumption (A4)
Let us consider the Markov process which is maybe the simplest one.That is, we consider that D = {1, 2} and X t jumps from 2 to 1 at rate a and cannot jump from 1 to 2. Similarly, it dies from 1 at rate b and cannot die from 1.It is a pure death process whose transition are represented in Figure 4. Its study is trivial and based on the spectral property of a 2 × 2 matrix.We have three cases : • b < a : We have exponential convergence to the unique QSD δ 2 .
• b = a : We have a slow convergence to the unique QSD δ 2 .• b > a : We have two QSD: and, furthermore, the convergence depends on the initial condition : if µ 0 ({2}) = η({2}) > 0 then the conditional law µ t converges to a b δ 1 + b−a b δ 2 .We see in the third case that δ 2 M t 1 ∼ e −at and δ 1 M t 1 ∼ e −bt are not of the same order and consequently (A4) fails and then the conclusion of Theorem 2 doesn't hold (which is what we naturally and simply observe here).
This process is trivial to study but it is an interesting example to understand the algorithms which approximate the quasi-stationary distribution.For instance the limit in large number of particles and in time of the well used Fleming-Viot algorithm do not commute for this example.Some question on its approximation remains open see [5].Commutation of limits for QSD approximation, as detailed in Section 5 below is then of prime importance.
A less simple and more interesting example is given by the house-of-card model studied in [19] and references therein.This process belongs to [0, 1] and is redrawn uniformly at random over this interval at rate 1.It is soft-killed at rate x → cx q , for some c, q ≥ 0. We can show the following : • If q > 1 − 1/c then we have exponential convergence to a unique explicit QSD (which has a bounded Lebesgue density).This in particular always verified when q ≥ 1 or c < 1.The convergence can be established in various distances with explicit rates.• If c = q+1 then we have two QSD : one, denoted by ν, has a Lebesgue density proportional to x → x −q and the second is δ 0 .When µ 0 ({0}) = 0, then we have a (explicit) slow convergence of the conditional laws µ t to the QSD ν (which has a Lebesgue density proportional to x → x −q ).When µ 0 ({0}) > 0 then µ t → δ 0 .Moreover, if q < 1/2 then η t e θt → c 0 ν, for some c 0 > 0 although if q ≥ 1/2 then η t e θt → 0. • If c < q + 1 then the conditional law converges, in mean, to the unique degenerate QSD: For details, additional results and references, see [19].Note that, in this example, Assumption (10) fails.

Lyapunov function and Assumption (A1) and (A2)
The Lyapunov assumptions require to find a kind of sub-solution and super-solution to the eigenvalue problem where the associated sub-eigenvalue β is strictly larger than the super-eigenvalue α.Let us show here a well-studied example of Markov processes where we easily find two Lyapunov functions in such a way that the condition α = β becomes a critical situation.
Let us consider the birth and death process on N * generated by for every n ≥ 2 and bounded function f , with b > d, and with some b 1 , d 1 > 0. It is known from [60] that e θt η t converge to a unique QSD if and only if This result can easily be recovered (and extended) by using Theorem 2 with V and ψ of the form q → q n with appropriate q > 0. See [3] for details.Note that the condition ( 10) is trivially satisfied.
3. The Langevin process (Round 1) In this section, we shall consider the following Langevin dynamics: Its infinitesimal generator is given for (x, v) ∈ R d × R d by: with formal adjoint L * in L 2 (dxdv): The metastable states often correspond to the basins of attraction of a potential V when the force F in (3) is conservative, i.e.F = −∇V .Therefore, we shall consider metastable states written as cylinders in the position and momenta coordinates as follows: where O is an open C 2 bounded connected set of R d .

Quasi-stationary distribution: existence and long-time convergence
The justification of the Parallel Replica method mentioned in Section 1.2 for the Langevin dynamics (13) in D first requires the existence of a QSD on D for the associated process (X t ) t≥0 killed outside of D. Such existence can be obtained through the application of the Krein-Rutman theorem to the semigroup of (X t ) t≥0 , which shall provide a left eigenvector of this semigroup corresponding to the QSD density.To apply the Krein-Rutman theorem, this semigroup is required to be compact in a Banach space.To prove the compactness the idea of the proof is to first show that the semigroup is an integral Hilbert-Schmidt operator in L 2 (D) hence compact in L 2 (D) using appropriate Gaussian upper-bounds for its transition density obtained in [49,Theorem 2.19].Second, this compactness is propagated across the Banach space C b (D) of continuous and bounded functions on D using the continuity of the semigroup, leading to the following result in [47, Theorem 2.9].
Theorem 3 (Compactness of the killed semigroup).For any t ≥ 0, p ∈ [1, +∞] and f ∈ L p (D), the quantity (2) For any t > 0, the operator M t is compact from L p (D) to L p (D), and from The application of the Krein-Rutman theorem to the operator M t for t > 0 on the Banach space C b (D) then yields the existence of a QSD on D of (X t ) t≥0 , with a density belonging to C b (D).The uniqueness of this QSD follows from the positivity of the transition density of the semigroup using a Harnack inequality provided in [49].In addition, one can also obtain a spectral interpretation of this QSD using Itô's formula, see the proofs of [47, Theorems 2.14, 2.17] for more details.
Theorem 4 (Existence and uniqueness of a QSD).There exists a unique QSD µ on D for the killed process (X t ) t≥0 .In addition, there exists a smooth positive function ψ such that where n(x) is the unitary outward normal vector at x ∈ ∂O.
In addition, the compactness of the killed semigroup allows to obtain a spectral decomposition of the killed semigroup.This gives in particular the long-time asymptotics of the law of the dynamics conditioned to remain in D, see [47,Theorem 2.22].
Theorem 5 (Convergence to the QSD in total variation).There exists α > 0 such that for any probability measure ν on D, there exists C ν > 0 such that for all t ≥ 0, where µ t is the law of (X t ) t≥0 conditioned to survival defined in (2) satisfying X 0 ∼ ν, and ∥ • ∥ T V denotes the total-variation norm on the space of bounded signed measures on R 2d .
The QSD on D can thus be seen as the longtime limit of the dynamics conditioned to stay inside D. This proposition is useful to understand what is a metastable state.A metastable state is a state such that the exit event from this state happens on a timescale larger than the local equilibration time, namely the time to observe the convergence to the QSD in (18).
Furthermore, if the Langevin dynamics ( 13) is initially distributed according to the QSD in D, then we can explicitly provide the law of the first exit event from D, which is given by the couple of the first exit time from D and its associated exit point.As mentioned in the introduction, recall that the first exit time from a domain when starting from a QSD is well known to be distributed according to an exponential law; this is a general property for Markov processes.By contrast, the law of the first exit point is strongly dependent on the dynamics.In the case of the Langevin dynamics (13), this is the first result to compute the law of the first exit point starting from the QSD.
Proposition 1.Let us assume that X 0 is distributed according to the QSD µ in D. Then the law of the couple (τ, X τ ) is fully characterized by the following properties: • τ follows the exponential law of parameter θ where θ is defined in Theorem 4; • τ is independent of X τ ; • The law of X τ is given by: for any bounded function f : where n(x) is the unitary outward normal vector at x ∈ ∂O and σ ∂O denotes the Lebesgue measure on the surface ∂O.
As a result, the existence of a QSD and the long-time convergence to the QSD of the Langevin dynamics trapped in D ensure that the Parallel Replica method also applies for the Langevin process trapped in D following the reasoning made in [45, Proposition 5].

Quasi-stationary distribution: overdamped limit
In this section we shall assume that γ > 0 and σ = 2γβ −1 for some β > 0 in (13).The aim of this section is to describe the behavior of the Langevin dynamics (13) when the friction coefficient γ goes to infinity.This will allow us in particular to explicit the overdamped limit of the QSD.
It is well known in the literature that the law of the rescaled position of the Langevin dynamics converges to the law of the overdamped Langevin dynamics (4), i.e. for all T > 0: However, this result does not provide an overdamped limit for the couple of the position and velocity vectors.Such a limit was provided recently in [58] where the following limit is obtained: where Z ∼ N d (0, β −1 I d ) is a Gaussian vector independent of the process (x t ) t∈[0,T ] .The proof of (20) consists in perturbing the Brownian noise of the Langevin dynamics ( 13) by a vanishing term when γ goes to infinity such that the position and velocity become independent.Then it remains to consider the overdamped limit of the position and velocity coordinates separately, which is much easier.In addition, since the perturbation converges to zero when γ goes to infinity, the overdamped limit of the perturbed process is the same as the overdamped limit of the Langevin dynamics.Following this result one can identify explicitly the overdamped limit of the QSD for the Langevin dynamics.
In particular, if one considers the marginal in position of the Langevin QSD then it converges weakly to the QSD of the overdamped dynamics on O as a consequence of the above proposition.

Langevin processes (round 2)
In this section, we present the main results of [36,37].These results first aim at giving necessary conditions for existence and uniqueness of a quasi-stationary distribution µ D of a strongly Feller Markov process (X t , t ≥ 0) killed when it exits a subdomain D, see Section 4.1.These results are then applied to Langevin processes with continuous or singular potentials, see Sections 4.2 and 4.3.

Notations and assumptions
Let (X t , t ≥ 0) be a time homogeneous Markov process valued in a Polish space S, with càdlàg paths and satisfying the strong Markov property, which is defined on the filtered probability space (Ω, F, (F t ) t≥0 , (P x ) x∈S ) where P x (X 0 = x) = 1 for all x ∈ S. Its transition probability semigroup is denoted by (P t , t ≥ 0).
Let B(S) be the Borel σ-algebra of S, bB(S) the space of all bounded and Borel measurable functions f on S (its norm will be denoted by f ∈ bB(S) → ∥f ∥ bB(S) = sup x∈S |f (x)|).The space D([0, T ], S) of S-valued càdlàg paths defined on [0, T ] is equipped with the Skorokhod topology.
We suppose that (C1) (strong Feller property) There exists t 0 > 0 such that for each t ≥ t 0 , P t is strong Feller, i.e.P t f is continuous on S for any f ∈ bB(S).(C2) (trajectory Feller property) For every T > 0, where D c = S\D is the complement of D. The transition semigroup of the killed process (X t , 0 ≤ t < τ ) is: Assume the following on the killed process (X t , 0 ≤ t < τ ): (C4) For t ≥ 0, M t is Feller, i.e. if f is bounded and continuous over D, so is M t f .(C5) There exists t 0 > 0 such that for all t ≥ t 0 , x ∈ D, and non-empty open subsets O of D, (we can assume this t 0 > 0 is the same as the one in (C1)), and there is some x 0 ∈ D such that P x0 (τ < +∞) > 0.

The general result for existence and uniqueness of the quasi-stationary distribution
In this section, we provide a slightly less general result than [36,Theorem 2.2].Indeed, Theorem 6 below is actually still valid when replacing (C4) by the less stringent assumption that for t ≥ 0, M t is weakly Feller.Theorem 6. Assume that (C1)→(C5) hold.Take any p ∈]1, +∞[.Then, there exists only one quasi-stationary distribution µ D of the process (X t , t ≥ 0) in D satisfying In addition, there are some constants δ > 0 and C ≥ 1 such that for any initial distribution ν on D with ν(W 1/p ) < +∞, where φ : D → R is continuous, φ/W 1/p is bounded, and φ is positive within D.
Let us give some remarks on this theorem.If W is bounded over D, the quasi-stationary distribution inside D is unique.When S is moreover compact, Theorem 6 is valid without assuming (C3) and with W = 1 there (i.e.there exists a unique quasi-stationary distribution in M 1 (D) and ( 23) holds for any ν ∈ M 1 (D)).On the contrary, when W is unbounded, Theorem 6 states that there is a unique QSD which integrate W 1/p for all p > 1 (in the statement the QSD µ D a priori depends on p, but then, clearly, it doesn't, since µ D (W 1/p ) < ∞ for some p > 1 implies µ D (W 1/q ) < ∞ for all q ⩾ p), but it doesn't prevent the existence of other QSD ν with ν(W 1/p ) = +∞ for all p > 1.This is for instance the case for the Ornstein-Uhlenbeck process killed at 0 (see [50]), to which Theorem 6 applies with W(x) = e αx 2 with α < 1/2.
The intuition behind the assumptions (C1)→(C5) is the following.It is now well known that a Lyapunov condition provides a spectral gap for a non killed semigroup in some weighted Banach spaces [29,52].Using results of [64] characterizing the essential spectral radius of a positive operator (combined with the regularity conditions (C1) and (C2) on the non killed process (X t , t ≥ 0)), the Lyapunov type condition (C3) is used to get a spectral gap for M t in b W 1/p B(D) (the space of measurable functions f over D such that f /W 1/p is bounded).It is also known that an irreducibility condition and a regularity condition on a non killed semigroup ensure the uniqueness of the invariant Notice that when (25) holds, Assumption (Ac1) is not satisfied but (Ac0) is satisfied.When V , c, and Σ satisfy respectively (Av1), (Ac0), and (AΣ), there is a unique weak solution to (24) by [63,Lemma 1.1], which is thus a (strong) Markov process.In the following we always assume that (AΣ) holds.
As explained in the introduction above (see also [27]), the process (x t , t ≥ 0) is stuck during long periods of time in neighborhoods of the local minima of V .This is due to energetic barriers this process has to cross to visit other regions of R d .Such neighborhoods are called metastable regions and are thus of the form where O ⊂ R d .We are interested in the existence of quasi-stationary distributions for the processes (24) in domain D of this form.We assume that O is a Some examples of functions V and c satisfying (Av1), (Av2), (Ac1), and (Ac2) are given in [63,Remark 3.2], see also [36,Remark 6.3].The Hamiltonian of the process ( 24) is, for Assume that (AΣ), (Av1), (Av2), (Ac1), and (Ac2) hold.Let us introduce for (x, v) ∈ R 2d , the modified Hamiltonian [63, Eq. (3. 3)], where G, U are as (Av2) and (Ac2), a > 0, b > 0, and w : R d → R is a compactly supported C 2 function.Define, for all x, v ∈ R d : The following asymptotic upper bound on W 1 gives a way to check if a probability measure on ν on D ⊂ R 2d satisfies ν(W 1/p ) < +∞ (see the proof of [36, Lemma 2.6]).
When O is bounded, V is C ∞ on O, and both Σ and c are constant, the very recent work [47] (see also [49,58]) establishes the uniqueness of the quasi-stationary distribution in the whole space of measure over D. These results were presented in Section 3 where the authors also considered non gradient vector fields (see indeed (13)).Their approach is different from the one adopted in [36,37]: it is based on a Feynman-Kac type formula for the killed semigroup, Harnack inequalities, and Gaussian upper bounds.We finally also refer to [4] for the recent study of Hypoelliptic diffusions killed at the boundary of a bounded subdomain of R m with non characteristic boundary.4.2.2.The case when (25) and (26) are satisfied When ( 25) and ( 26) hold, and ℓ 0 < n 0 −2, we are able to construct a bounded Lyapunov function for the process (24) (i.e. a bounded Lyapunov function W over R 2d satisfying (C3)).This implies the following theorem.
Then, Theorem 6 is valid with a bounded Lyapunov function.In particular: there exists a unique quasi-stationary distribution in D = O × R d for the process (24) in the whole space of probability measures M 1 (D) on D, and in addition, (23) holds for all ν ∈ M 1 (D).

Langevin processes with singular potentials
In this section we consider the existence and uniqueness of a quasi-stationary distribution for the process (24) when the potential V is singular.For ease of exposition, we will just focus here on Lennard-Jones and Coulomb interactions.However much more general singular potentials are considered in [37].
For d ≥ 1, consider a system of N ≥ 2 particles in R d which cannot collide and let denote respectively the positions of the N particles and their velocities, at time t ≥ 0. The natural space to consider the time evolution of the positions (x t , t ≥ 0) and of the velocities (v t , t ≥ 0) of the N particles is thus Notice that in both cases, O is open, path connected, and unbounded.In molecular dynamics, the interatomic potential of the system of N particles is typically of the form, for x = (x 1 , . . ., x N ) ∈ O, where V c : R d → R is the confining potential of the system and V I : Ω → R (where, if d = 1, Ω = {y < 0}, and if d ≥ 2, Ω = {y ̸ = 0}) is the potential energy modeling the interaction between two particles, the latter becoming infinite when (and only when) y ∈ ∂Ω = {0} (which prevents from collisions).Let (Ω, F, (F t ) t≥0 , P) be a filtered probability space.We assume that the evolution of the positions (x t , t ≥ 0) and the velocities (v t , t ≥ 0) of the N particles on S is described by the following hypoelliptic stochastic differential equation where c We set X t = (x t , v t ) for t ≥ 0. Throughout this section, we assume that Σ satisfies (AΣ) (with N d instead of d there) and that c is such that: In the two next sections, we consider existence and uniqueness of a quasi-stationary distribution for the process (32) when V I is the Lennard-Jones potential or the Coulomb potential.We recall that the Lennard-Jones potential is defined by, for y ∈ Ω, where b, c > 0. When d ≥ 3, the Coulomb potential is the function defined on Ω by: where e > 0. In this section, the interaction potential V I is either V co or V LJ .The Coulomb potential when d = 2 is treated in the next section.We assume that V c is C 2 over R d and to simplify we assume that for some r > 0, where A > 0 and α 1 > 1.We refer to [37] for more general confining potentials.With such confining potentials and interaction potentials, by [37, Proposition 2.3], for any x ∈ S, there exists a unique pathwise solution (X t = (x t , v t ), t ≥ 0) of ( 32) with X 0 = x, which is moreover non-explosive and remains in S for all t ≥ 0.
In other words, Theorem 9 states that on D, there exists a unique quasi-stationary distribution for the process (32) in the space M p (∀p > 1).In addition, Equation (23) holds for all ν ∈ M p .
We refer to [37,Proposition 2.10] for the explicit construction of W, which is inspired by the previous works [63] and [39].Notice that O is not necessarily bounded in Theorem 9, and its closure may contain singularities of V , namely some subset of ∂O.In this section, we consider the interaction potential V I defined by for all y ∈ R d : V When d = 2, V I is the Coulomb potential, and when d = 1, V I corresponds to a log singularity pairwise potential.We also assume that the confining potential ) and for some r > 0, it holds: where A > 0 and α 2 ≥ 0 (we refer to [37] for more general confining potentials).
With such an interaction potential and confining potentials, by [37, Proposition 3.1], for any x ∈ S (see (30)), there exists a unique pathwise solution (X t = (x t , v t ), t ≥ 0) of (32) with X 0 = x, which is moreover non-explosive and remains in S (see (30)) for all t ≥ 0. Theorem 10.Let d ∈ {1, 2}.Assume that (Ac) and (AΣ) are satisfied.Assume also that V I is given by (35) and that V c is a C 2 potential such that (36) Then, Theorem 6 is valid for the process (32) on D with a continuous and unbounded Lyapunov functional W : S → [1, +∞) which satisfies, for some m > 0, W ≤ exp mH η2 on S.
We refer to [37,Proposition 3.3] for the explicit construction of W, which is inspired by [51].Let us emphasize that there is no restriction in Theorem 10 on η 2 (i.e. one can choose any η 2 in (0, 1]).We mention that in [37], we also deduce that large deviation principles hold for the non killed process (i.e. the process (X t , t ≥ 0) on S) with Lennard-Jones and Coulomb potential interactions, as well as the exponential convergence of its law towards its invariant measure.The existence and uniqueness of a quasi-stationary distribution for elliptic diffusions with Lennard-Jones type interactions is also studied there.

Numerical approximation with particles
This section is based on [40].The goal is to introduce a numerical scheme whose aim is to sample the QSD of a killed diffusion, and to prove quantitative long-time convergence rates for this algorithm.

The problem
Given a Markov process X killed at time τ , we introduce a numerical scheme in order to sample its quasi-stationary distribution ν * , or the law of X t conditioned on survival.With great generality, one may introduce a particle approximation of the process conditioned on survival, called a Fleming-Viot process in most of the literature.It consists in N ∈ N independent processes (X i ), having the same dynamic as X, but when one of them dies, instead of staying at the cemetery point, it chooses uniformly at random another one of the processes, and branches onto it.This gives rise to a mean-field interacting particle system, with Markov semi-group (P N,t ) t⩾0 .In the particular case of a softly-killed elliptic SDE on the d-dimensional torus: where b ∈ C 1 (T d ) and B is a d-dimensional Brownian motion, we propose a numerical scheme of this Fleming-Viot process.Given a random variable E ∈ R + with law e −t dt, we define the soft killing as: where λ : T d → R + is a L λ -Lipschitz death rate, for some L λ > 0. As before, write: the law of the process conditioned on survival.First, we introduce a discretized version of the killed diffusion.We define a Markov kernel K γ on T d for γ small enough as follow: fix x ∈ T d , then We define a Markov chain (X n ) by X n+1 ∼ K γ (X n , dx).Given a sequence of independent and uniform on [0, 1] random variable (U n ), we define the death time of the chain (X n ) as: where p(x) = 1 − e −γλ(x) .Write η n = Law(X n |τ > n).The discretized killing as been defined as such in order to have for all t ⩾ 0: as γ → 0. This Markov chain also admits a QSD ν γ .Now, we want to define two processes: a Markov process (Y n ) whose law is η n , defined for all time, and a particle approximation of this Markov chain (Y n ).To this end, we define first, given some probability measure µ ∈ M 1 (T d ), a Markov kernel Q µ such that Q µ (x, •) corresponds to the law of some random variable having done a step of discretized diffusion, and upon death, it does a step of discretized diffusion, conditioned on survival, with initial condition µ.More precisely, we construct a random variable of law Q µ (x, •) as follow: (1) Draw X 0 ∼ K γ (x, •) and U 0 ∼ U([0, 1]).
(2) If U 0 ⩾ p(X 0 ), set X = X 0 in T d (in that case, we say the particle has moved from x to X 0 without dying).(3) If U 0 < p(X 0 ) then set i = 1 and, while X is not defined, do: (a) Draw a new X ′ i distributed according to µ, a new X i ∼ N (X ′ i + γb(X ′ i ), γI d ) and a new U i ∼ U([0, 1]).(b) If U i ⩾ p(X i ), set X = X i in T d (in that case, we say the particle has died, resurrected at X ′ i , moved to X i and survived).(c) If U i < p(X i ), set i ← i + 1 (in that case, we say the particle has died, resurrected at X ′ i , moved to X i and died again) and go back to step (a).We introduce the following non-linear process, solution to the problem: This non-linear Markov chain does at each time n a step of discretized diffusion conditioned on survival.As said, the reason to introduce this Markov chain is the following property: Now we can define the numerical scheme of the Fleming-Viot process.For x = (x 1 , . . ., x N ), we write: the empirical measure and: We then define recursively the process by X n+1 ∼ R γ,N (X n , •).This can be understood as follow: each particles is an Euler-Maruyama scheme of the diffusion, and if one of them dies at time n, then it does a step of K γ , starting from π(X n−1 ), conditioned on survival.As N goes to infinity, if the initial condition condition of this particle system is N iid random variable of law η 0 , then π(X n ) is expected to be close to the common law of the X i 's.This can then be shown to be close to η n .This is the propagation of chaos phenomenon.

Main results
For this Markov chain, we proved two theorem.The first one is about the long-time behavior of the process: Theorem 11.There exists c 0 , γ 0 > 0 that depends only on the drift b and the dimension d such that, if λ is Lipschitz with a constant L λ such that then there exist a distance ρ equivalent to | • | on T d , such that ρ N (x, y) = ρ(x i , y i ) is a distance on T dN , and κ > 0, such that for all γ ∈ (0, γ 0 ] N ∈ N, and all x, y ∈ T dN , there exists a coupling (X, Y) of δ x R N,γ and δ y R N,γ such that: This yields the existence of a stationary measure µ ∞,N,γ , and exponential convergence of the Markov chain towards this probability measure.
To construct the coupling of theorem 11, we first construct a coupling of Q µ (x, •) and Q ν (x, •), for any x, y ∈ T, µ, ν ∈ P(T d ).This coupling is done as follow: thanks to [30], we are able couple the step of discretized diffusion to get a coupling (X ′ , Y ′ ) of δ x K γ and δ y K γ .The probability that one of them dies and not the other is less then L λ E(|X ′ − Y ′ |).If they die together, then we construct again thanks to [30] a sequence of coupling of µK γ and νK γ .
The second theorem is about the convergence of the empirical measure of the process towards the quasi-stationary distribution of the process (37) killed at time τ , in the limit γ → 0 and n, N → ∞.Theorem 12.Under the hypothesis of theorem 11, there exists C > 0 such that for all N ∈ N, γ ∈ (0, γ 0 ], t ⩾ 0 and µ 0 ∈ P(T dN ), if (X k ) k∈N is a Markov chain with initial distribution µ 0 and transition kernel R N,γ , where The second theorem is a consequence of the first one.The rate of convergence in γ, N and t would be the same in the case of N independent diffusion, hence are optimal.The rate of convergence in N comes from [33].
To prove propagation of chaos, we construct a coupling of the discretized Fleming-Viot process, with a system of N independent non-linear process defined in equation (40).This coupling is done is a similare fashion as for the previous coupling of R γ,N (x, •) and R γ,N (x, •).For the limit γ → 0, the coupling of the process is made with a continuous time Fleming-Viot process.Actually, thanks to those couplings, we are able to prove similar convergence when only one or two parameters among t, γ, or N are sent to their limit.We can summarize our results thanks to the diagram of figure 5.

Figure 1 .
Figure 1.A sample path of the chemostat model.Here the blue curve corresponds to the discrete component (N t ) and the green curve corresponds to the continuous component (S t ).

Figure 3 .
Figure 3.A sample path of the process generated by (12).

Figure 4 .
Figure 4.The transition diagram of the Markov chain described in Subsection 2.2.2 ) is well-defined.Besides,(1) the family of operators (M t ) t≥0 is a semigroup on L p (D) and on C b (D).

where 1
continuous from S to the space M 1 (D([0, T ], S)) of probability measures on D([0, T ], S), equipped with the weak convergence topology.(C3) (Lyapunov function condition) There exist a continuous function W : S → [1, +∞[, belonging to the extended domain D e (L) of the generator of (X t , t ≥ 0), two sequences of positive constants (r n ) and (b n ) where r n → +∞, and an increasing sequence of compact subsets (K n ) of S, such that −LW(x) ≥ r n W(x) − b n 1 Kn (x), q.e, Kn is the indicator function of K n .Now let D be a non-empty open domain of S, different from S. Consider the first exit time of D τ := inf{t ≥ 0, X t ∈ D c }

Figure 5 .
Figure 5. Summary of the results