A M G Cox

E Horton
email: emma.horton@inria.fr

D Villemonais
email: denis.villemonais@univ-lorraine.fr

Binary branching processes with Moran type interactions

Keywords: interacting particle systems, branching processes, many-to-one, Markov processes, birth-and-death process, Moran model 82C22, 82C80, 65C05, 60J25, 92D25, 60J80

 and population size dynamics.

Introduction

Branching processes naturally arise as pertinent models in population size dynamics [START_REF] Jagers | General branching processes as Markov fields[END_REF][START_REF] Jagers | Branching processes as population dynamics[END_REF][START_REF] Lambert | The branching process with logistic growth[END_REF], neutron transport [START_REF] Cox | Multi-species neutron transport equation[END_REF], genetic dynamics [START_REF] Marshall | A branching process for the early spread of a transposable element in a diploid population[END_REF], growth-fragmentation processes [START_REF] Bertoin | Markovian growth-fragmentation processes[END_REF][START_REF] Bertoin | On a Feynman-Kac approach to growth-fragmentation semigroups and their asymptotic behaviors[END_REF] and cell proliferation kinetics [START_REF] Yanev | Branching processes in cell proliferation kinetics[END_REF], and as theoretical objects in their own right [START_REF] Jagers | General branching processes as Markov fields[END_REF][START_REF] Dawson | Measure-valued Markov processes[END_REF][START_REF] Ikeda | Branching Markov processes I[END_REF][START_REF] Ikeda | Branching Markov processes II[END_REF][START_REF] Ikeda | Branching Markov processes III[END_REF]. These models are characterised by the independence of the branching and killing events in the system, which leads to a multiplicative behaviour, as well as interest in the scaling behaviour and the characterisation of the system at large times. On the other hand, processes with Moran type interactions are natural models for finite populations with either variety-increasing or variety-reducing effects such as genetic drift, genetic mutations and natural selection. First introduced by Moran [START_REF] Moran | Random processes in genetics[END_REF], the Moran model describes the evolution of N genes such that at an exponential rate, two particles are chosen uniformly at random from the population, one of which is killed and the other splits in two. Thus, the independence between particles is lost. We refer the reader to [START_REF] Etheridge | Some Mathematical Models from Population Genetics: École D' Été de[END_REF] and references therein for an overview of this model and of its extensions. This type of resampling has since been employed in a range of particle system models to numerically solve Feynman-Kac models, see [START_REF] Cloez | Uniform in time propagation of chaos for a Moran model[END_REF][START_REF] Del Moral | A Moran particle system approximation of Feynman-Kac formulae[END_REF][START_REF] Del Moral | Branching and interacting particle systems approximations of feynman-kac formulae with applications to non-linear filtering[END_REF][START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF] and references therein.

The model we consider in this paper provides a combination of these two types of dynamics: (natural) branching and killing, as well as Moran type interactions. More precisely, when the system is initiated from N particles, each particle evolves according to an independent copy of a given Markov process, X, until either a (binary) branching or killing event happens. Here, binary refers to the fact that the particle is replaced by exactly two independent copies of itself. If such a branching event occurs, with a probability that may depend on the configuration of the whole system, another particle is removed from the system a selection mechanism. Similarly, if a killing event occurs, with a probability which may also depend on the configuration of the whole system, another particle is duplicated via a resampling mechanism. We refer to this model as the binary branching model with Moran interactions, or BBMMI for short.

Our main result (cf. Theorem 6) gives an explicit relation between the average of the empirical distribution of the particle system and the average of the underlying Markov process X. Indeed, letting m T denote the empirical distribution of the particle system at time T , we show that for any T ≥ 0,

E Π A T Π B T m T f = E m 0 f (X t ) exp t 0 b(X s ) -κ(X s ) ds 1 t<τ ∂ , (1)
where b is the branching rate, κ is the soft killing rate, τ ∂ is the absorption time, Π A T and Π B T are weights that compensate for the resampling and selection events that occur up to time T , which we state explicitly in Theorem 6, and E • is the expectation with respect to the law of X.

As we will later discuss, by choosing particular model parameters, we recover both the classical many-to-one formula [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF] and the unbiased estimator proved in [START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF] for fixed size Moran type genetic algorithms. The second part of Theorem 6 then provides a precise bound on the L 2 distance between the normalised semigroup on the right-hand side of (1) and the normalised occupation measure of the particle system. In particular, it states that this distance converges to zero as the number of initial particles tends to infinity.

In the third section of this article, we illustrate several applications of this model by focussing on a particular selection/resampling mechanism in our model. More precisely, we consider a binary branching process whose population size is constrained to remain in {N min , . . . , N max }, where 2 ≤ N min ≤ N max < +∞. In order to constrain the size of the process, when the size of the population reaches N max (resp. N min ) and a natural branching (resp. killing) event occurs, we set the probability of selection (resp. resampling) to be 1. We often refer to this model as the N min -N max model. As will be clear in the next section, the choice N min = 0 and N max = +∞ is also allowed, so that the collection of models we consider ranges from constant population size models with neither branching nor killing (corresponding to the choice N min = N max ) to branching population models (when N min = 0 and N max = +∞) with or without a selection/resampling mechanism.

For comparable results concerning the resampling algorithm, we refer the reader to [START_REF] Ferrari | Quasi stationary distributions and Fleming-Viot processes in countable spaces[END_REF][START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF][START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF][START_REF] Oçafrain | Convergence of a non-failable mean-field particle system[END_REF] and to [START_REF] Delyon | A central limit theorem for Fleming-Viot particle systems with hard killing[END_REF][START_REF] Del Moral | Asymptotic results for genetic algorithms with applications to nonlinear estimation[END_REF][START_REF] Lelièvre | Central limit theorem for stationary Fleming-Viot particle systems in finite spaces[END_REF] for associated central limit theorems. Whether a central limit theorem can be proved in our setting remains open, but could be approached by carefully studying the square increments of the martingale decomposition used in the proof of Theorem 6.

Our model and results are also reminiscent of the genetic algorithms introduced by Del Moral (see [START_REF] Moral | Feynman-Kac formulae. Probability and its Applications[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF] and references therein), where a particle system with killing is constrained to remain at a constant size via a resampling mechanism. We also refer the reader to [START_REF] Burdzy | Configurational transition in a Fleming-Viot-type model and probabilistic interpretation of Laplacian eigenfunctions[END_REF][START_REF] Burdzy | A Fleming-Viot particle representation of the Dirichlet Laplacian[END_REF][START_REF] Grigorescu | Hydrodynamic limit for a Fleming-Viot type system[END_REF][START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF][START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF][START_REF] Cloez | Uniform in time propagation of chaos for a Moran model[END_REF] and references therein.

Finally, we note that our model fits into the more general class of controlled branching processes introduced by Sevastyanov and Zubkov in [START_REF] Sevastyanov | Controlled branching processes[END_REF], where the number of reproductive individuals in each generation depends on the size of the previous generation via a control function. This work was later extended by Yanev [START_REF] Yanev | Conditions for degeneracy of ϕ-branching processes with random ϕ[END_REF] to allow for random control functions. We refer the reader to [START_REF] Sevastyanov | Controlled branching processes[END_REF][START_REF] González | On the class of controlled branching processes with random control functions[END_REF][START_REF] Velasco | Controlled branching processes[END_REF][START_REF] Yanev | Critical controlled branching processes and their relatives[END_REF] for discrete time versions of this process and to [START_REF] Garcia | Continuous-time controlled branching processes[END_REF] for the continuous time version. In these articles, the authors study the convergence of the survival probability in different regimes, as well as the expected population size. We also note that the control function can be seen as a way to model immigration and emigration of particles, [START_REF] Yanev | Critical controlled branching processes and their relatives[END_REF]. We refer the reader to [START_REF] Vatutin | A critical Galton-Watson branching process with emigration[END_REF][START_REF] Kawazu | Branching processes with immigration and related limit theorems[END_REF][START_REF] Olofsson | General branching processes with immigration[END_REF][START_REF] Li | Branching processes with immigration and related topics[END_REF][START_REF] Lambert | Population dynamics and random genealogies[END_REF][START_REF] Barczy | A note on asymptotic behavior of critical Galton-Watson processes with immigration[END_REF][START_REF] Bansaye | Spine for interacting populations and sampling[END_REF] for results regarding these latter processes.

The rest of the paper is set out as follows. In the next section, we introduce some notation and assumptions, and formally describe the particle system introduced above. We will then state our main result in Theorem 6 and discuss several implications.

In section 3, we focus on the N min -N max model and consider a variety of applications in order to illustrate the scope of our results, as well as some of the possible difficulties and extensions. In particular, in Section 3.1, we describe and study a stochastic population model with constrained size, based on the N min -N max model. We give a sufficient condition ensuring that the number of resampling events does not explode in finite time when the killing rate of X is unbounded. We also relate the long-time convergence of the empirical distribution of the population to the existence of a quasi-stationary distribution for an auxiliary sub-Markov process related to X and show the time uniform convergence of the renormalised empirical distribution of the particle system toward its conditional distribution, at a polynomial speed when N max /N min goes to infinity. We prove these results in the particular setting of this model, but our methods could be extended to a more general setting.

In Section 3.2, we deal with the particular case where the underlying Markov process is a piecewise deterministic Markov process (PDMP). In this case, it is possible for two particles to be killed at the same time, meaning that our main result is not directly applicable without further work. Working in the particular setting of the neutron transport equation (NTE), see [START_REF] Horton | Stochastic methods for the neutron transport equation I: Linear semigroup asymptotics[END_REF][START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF], we overcome this difficulty using the notion of h-transform introduced in [START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF] for the NTE. We show that particles whose dynamics are given by an appropriately h-transformed process fulfil the assumptions of our main theorem, and thus, we may use the N min -N max model with the h-transformed process in order to obtain estimates for quantities associated to the original process.

In Section 3.3, we focus on some numerical properties of the resampling/selection process. First, we make comparisons between the N min -N max model and the fixed size Moran type interacting particle system (IPS) for the approximation of quantities related to non-conservative semigroups (cf. [START_REF] Moral | Feynman-Kac formulae. Probability and its Applications[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF][START_REF] Burdzy | Configurational transition in a Fleming-Viot-type model and probabilistic interpretation of Laplacian eigenfunctions[END_REF][START_REF] Burdzy | A Fleming-Viot particle representation of the Dirichlet Laplacian[END_REF]), which corresponds to the particular case N min = N max with no branching. We observe that in the case of a birth-and-death process, X, with branching rate given by X t ∧ M for some M > 0 and state-dependent Markov transitions, the computational cost, measured as the number of resampling/selection events, is significantly lower when using the N min -N max algorithm. We also study the bias and standard deviation of an estimator for the normalised empirical stationary distribution for each of the N min -N max and fixed size Moran type IPS. Letting M go to infinity, we provide numerical evidence that the number of interactions, bias and standard deviation stabilise for the N min -N max process, whereas these quantities grow linearly for the fixed size Moran type IPS process. Based on these simulations, we conjecture that the results of Section 2.2 hold true in situations where the branching rate is unbounded and thus where the Moran type process is undefined. In the final part of Section 3.3, we introduce a filtering method to deal with some of the numerical discrepancies that occur when the variance of the number of resampling/selection events is high.

In Section 4 we provide the proof of Theorem 6. Finally, we end the paper with an appendix containing the formal construction of the particle system introduced in section 2.1.

Main results

Description of the model

Let (Ω, F, (X t ) t∈[0,+∞) ) be a continuous time progressively measurable Markov process with values in a measurable state space E ∪ ∂, where ∂ is an absorbing (measurable) set such that ∂ ∩ E = ∅. Denoting by τ ∂ = inf{t ≥ 0, X t ∈ ∂} its absorption time, it follows that

X t ∈ ∂, ∀t ≥ τ ∂ .
We denote by P x its law when initiated at x ∈ E and by E x the corresponding expectation operator. We extend, whenever necessary, any measurable function f : E → [0, +∞) by f ≡ 0 on ∂. We also assume that we are given two bounded functions b : E → R + and κ : E → R + , denoting.

We will shortly define a particle system, where particles move according to copies of X between interactions. Before doing so, we make the following assumption, necessary for the system to be well defined. It entails that two independent copies of X are absorbed simultaneously with probability 0. Assumption 1. For any x ∈ E and t ∈ [0, +∞), P x (τ ∂ = t) = 0 and P x (τ ∂ > t) > 0.

We also introduce the following notation. Denote by P f (N) the collection of finite subsets of N := {1, 2, . . . }. Then, for i ∈ N, we let

b i : (E ∪ ∂) P f (N) → [0, +∞), κ i : (E ∪ ∂) P f (N) → [0, +∞)
be a collection of bounded measurable functions. Also for i ∈ N, let

p i : (E ∪ ∂) P f (N) → [0, 1], q i : (E ∪ ∂) P f (N) → [0, 1]
be measurable functions. We will assume that, for any s ∈ P f (N), for any collection of points {x i : i ∈ s, x i ∈ E}, and any

i 0 ∈ s, b i 0 (x i , i ∈ s) -κ i 0 (x i , i ∈ s) = b(x i 0 ) -κ(x i 0 ) (2)
and that

p i 0 (x i , i ∈ s) = 0 whenever |s| = 1, (3)
where |s| denotes the number of elements in the set s.

We now proceed with an informal description of the dynamic of the particle system, which we call the binary branching model with Moran type interactions (BBMMI). The formal construction of the process is a non-trivial task, and is given in the appendix. To this end, fix N0 ≥ 2. We consider the particle system ( St , (X i t ) i∈ St ) t∈[0,+∞) , where the component St ∈ P f (N) is the set enumerating the particles in the system at time t, initiated with S0 = {1, 2, . . . , N0 } and X i 0 = x i ∈ E for all i ∈ S0 . The number of particles in the system at any time t is denoted by Nt = | St |.

Evolution of the BBMMI.

1. The particles X i , i ∈ S0 , evolve as independent copies of X, and we consider the following times:

τ b,i 1 := inf{t ≥ 0, t 0 b i (X i s , i ∈ S0 ) ds ≥ e i 1 }, and
τ κ,i 1 := inf{t ≥ 0, t 0 κ i (X i s , i ∈ S0 ) ds ≥ E i 1 } and τ ∂,i 1 := inf{t ≥ 0, X i t ∈ ∂},
where e i 1 , E i 1 , i = 1, . . . , N0 are exponential random variables with parameter 1, and are independent of each other and everything else.

2. Denoting by i 0 the index of the (unique) particle such that τ b,i 0

1 ∧ τ κ,i 0 1 ∧ τ ∂,i 0 1 = τ 1 , where τ 1 = min i∈ S0 τ b,i 1 ∧ τ κ,i 1 ∧ τ ∂,i 1 , we set St = S0 for all t < τ 1 and (a) if τ 1 = τ b,i 0 1 , then a branching event occurs; (b) if τ 1 = τ κ,i
1 , then a soft killing event occurs; (c) if τ 1 = τ ∂,i 0 1 , then a hard killing event occurs.

3. Then a resampling or selection event may occur, depending on the following situations:

killing: if a (hard or soft) killing event occurred at the preceding step, then we say that i 0 is killed at time τ 1 and

• with probability p i 0 (X i τ 1 , i ∈ S0 ), the particle i 0 is removed from the system and a resampling event occurs: one chooses j 0 uniformly in S 0 \ {i 0 } and sets

X max S 0 +1 τ 1 = X j 0 τ 1 and Sτ 1 := S0 \ {i 0 } ∪ {max S0 + 1};
observe that the number of particles in the system at time τ 1 is then Nτ 1 = N0 ; we say that j 0 is duplicated at time τ 1 . • with probability 1-p i 0 (X i τ 1 , i ∈ S0 ), the particle i 0 is removed from the system; the set of particles at time τ 1 is enumerated by Sτ 1 := S0 \ {i 0 } and the number of particles in the system is then Nτ 1 = N0 -1; branching: if a branching event occurred at the preceding step, then we say that i 0 has branched at time τ 1 and

• with probability q i 0 (X i τ 1 , i ∈ S0 ), a new particle is added to the system at position X i 0 τ 1 and a selection event occurs: one chooses j 0 at random uniformly in S 0 ∪ {max S 0 + 1} and removes the particle j 0 from the system:

X max S 0 +1 τ 1 = X i 0 τ 1 and Sτ 1 := {max S0 + 1} ∪ S0 \ {j 0 };
observe that the number of particles at time τ 1 is thus Nτ 1 = N0 ; we say that j 0 is removed at time τ 1 and that max S 0 + 1 is born at time τ 1 ; • with probability 1 -q i 0 (X i τ 1 , i ∈ S0 ), a new particle is added to the system at position X i 0 τ 1 :

X max S 0 +1 τ 1
= X i 0 τ 1 and Sτ 1 = S0 ∪ {max S0 }, and we say that max S 0 + 1 is born at time τ 1 .

After time τ 1 the system evolves as independent copies of X until the next killing/branching event, denoted by τ 2 , and at time τ 2 it may undergo a resampling/selection event as above. By iteration, we define the sequence

τ 0 := 0 < τ 1 < τ 2 < • • • < τ n < • • • .
(a) Resampling. Just before time 0.05, the orange particle is killed but no resampling occurs. Around time 0.15, the blue particle is killed and the green particle is resampled to produce the red particle.

(b) Selection. At time 0.02, the purple particle branches, producing the orange particle but no selection occurs. At the second dashed line, the dark blue particle branches to produce the red particle and the light blue particle is selected and removed from the system. We will also make use of the following assumption, which ensures that the process described above is well defined at any time t ≥ 0. Assumption 2. The sequence (τ n ) n∈N converges to +∞ almost surely.

Let us now make several remarks concerning the above particle system and potential parameter choices.

Remark 1. The enumeration of the particle system does not keep track of the genealogy of the particles in the system in order to simplify notation; this could be done using the Ulam tree as a set of indices.

Remark 2. The functions b i 0 and κ i 0 are the rate at which the particle i in the system branches and is killed respectively (in addition to hard killing events occurring when a particle hits the absorbing set ∂). The simplest situation is when b i 0 (x i , i ∈ s) = b(x i 0 ) and κ i 0 (x i , i ∈ s) = κ(x i 0 ), but other configurations are possible. We emphasise that the balance condition ( 2) is crucial to ensure the relation with the semigroup Q defined in the next section, and should not be considered a technical assumption.

Remark 3. The functions p i 0 (respectively q i 0 ) denote the probability that a killing event (respectively a branching event) triggers a branching (respectively a killing) among the particles in the system. We call such an event a resampling event (respectively a selection event). In the context of Moran type models, they correspond to selection by death events (respectively selection by birth events). Condition (3) prevents a resampling event from occurring when the population decreases to 0 particles. We emphasize that, since ∂ is a set, p i 0 may depend on the killing location.

Remark 4. The above dynamic allows one to constrain the size of the process to remain between two bounds N min = 1 and N max , with 0 ≤ N min ≤ N max . In order to do so, one simply chooses

p i 0 (x i , i ∈ s) = 1 {|s|=N min } and q i 0 (x i , i ∈ s) = 1 {|s|=Nmax} .
In what follows, we refer to this dynamic as the N min -N max process. Other natural choices include, for instance,

p i 0 (x i , i ∈ s) = 1 s + 1 1 |s|≥2 and q i 0 (x i , i ∈ s) = 1 - 1 s + 1 ,
i.e. resamplings become less probable when the number of particles increases and selection becomes stronger when the number of particles increases, respectively. Note also that, since X is only assumed to be Markov, it can carry the time as a component, which would allow N min and N max to vary with time (in the former situation) or the selection pressure to vary with time (in the later example), by simply letting p, q depend on time.

Remark 5. In our model, we consider general killing, which corresponds to a conservative Markov process X killed at a possibly unbounded rate. Finding general criteria ensuring that Assumptions 1 and 2 hold is usually involved (it is trivial when the killing rate is bounded and there is no hard killing). In Section 3.1 we provide a criterion based on a Foster-Lyapunov type assumption ensuring the non-explosion of the system when the killing rate κ is unbounded (and when there is no hard killing) and also show in Section 3.2 that Assumption 2 holds true for a specific PDMP model, obtained as the h transform of a Neutron Random Walk killed at the boundary of domain. There is also a substantial literature on these types of problems in the presence of hard killing for the more classical fixed size Moran type systems, and, in particular, the interested reader may find it useful to adapt the methods developed in [START_REF] Burdzy | A Fleming-Viot particle representation of the Dirichlet Laplacian[END_REF][START_REF] Grigorescu | Hydrodynamic limit for a Fleming-Viot type system[END_REF][START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF][START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF] to our setting.

In the rest of this article, we denote by σ n and ρ n the n th selection and resampling events, respectively. We let A t denote the total number of resampling events up to time t, B t denote the total number of selection events up to time t and we set C t = inf{n ≥ 0 : t < τ n+1 } to be the total number of events up to time t. We will also use the following quantities throughout,

Π A t := At n=1 Nρn -1 Nρn , Π B t := Bt i=1 Nσn + 1 Nσn . ( 4
)
2.2 L 2 bounds for the empirical distribution of the process For all bounded measurable functions f : E → R, all t ≥ 0 and all x ∈ E, we set

Q t f (x) = E x f (X t ) exp t 0 (b(X s ) -κ(X s )) ds 1 t<τ ∂ ,
where we recall that τ ∂ = inf{t ≥ 0 : X t ∈ ∂}. This defines a Feynman-Kac semigroup (Q t ) t≥0 , which is related to the binary branching model where particles move as copies of X that are killed at rate κ and branch at rate b resulting in the creation of two independent copies of the original particle. The relation between Q and this process is given by the well known many-to-one formula (see for instance [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF] and references therein, see also Remark 8). In this section, we consider the BBMMI defined above and study its relation with Q.

In what follows, we let m t (resp mt ) denote the occupation measure (resp. normalised occupation measure) of the particle system,

m t = i∈ St δ X i t , t ≥ 0, mt = 1 N0 m t . (5)
The following result first provides an equality between a penalised version of the occupation measure of the system at any time T and the semigroup Q at time T . The second part is dedicated to the control of the L 2 distance between renormalised versions of the semigroup Q and the occupation measure of the process.

Theorem 6. Under Assumptions 1 and 2, the BBMMI satisfies, for all time T ≥ 0 and all bounded measurable functions f : E → R, the following many-to-one formula

m 0 Q T f = E Π A T Π B T m T f . (6)
Moreover,

m 0 Q T f m 0 Q T 1 E - m T (f ) m T (1 E ) 1 m T =0 2 ≤ C exp(c b ∞ T ) f ∞ m 0 Q T 1 E / N0 1 N0 , (7)
where C, c are positive constants.

The proof of this result is given in Section 4 but let us first make some further comments on the model and the above result.

Remark 7. In Theorem 6, m 0 is assumed deterministic with total mass N0 . The result immediately extends to the case where m 0 is random by taking the expectation in ( 6) and (7) (in the latter, the right hand term may be infinite). Moreover, the expressions may also be written in terms of the normalised measure, m, in which case (6) becomes m0

Q T f = E Π A T Π B
T mT f and (7) can be written as

m0 Q T f m0 Q T 1 E - mT (f ) mT (1 E ) 1 mT =0 2 ≤ C exp(c b ∞ T ) f ∞ m0 Q T 1 E 1 N0 ,
making explicit the 'classical' dependence N -1/2 0 on the number of particles which appears on the right-hand side.

Remark 8. When p = q ≡ 0, the particle system has the dynamic of the classical binary branching process described just before the statement of the theorem. In this situation, we have Π A T = Π B T = 1 almost surely, and we thus recover the classical many-to-one formula.

Remark 9. Although we have remained in the setting of binary branching in this article, it should be possible to extend the definition of the BBMMI to allow for more general offspring distributions when a (natural) branching event occurs by modifying the selection mechanism. We intend to address this in future work.

Remark 10. The definition of our model allows a population with varying size by including branching, death and Moran type interactions, which are key features in the context of population dynamics. However, when restricted to the situation where p = q ≡ 1 (so that the size of the process remains constant equal to N0 ), Theorem 6 extends the results of [START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF] by allowing selection events. Similar results were also obtained in [START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF] (see also references therein and the recent paper [START_REF] Cloez | Uniform in time propagation of chaos for a Moran model[END_REF] for an adaptation to discrete state spaces, with more general interaction mechanisms), where the process is constrained to remain of constant size N0 and with uniformly bounded killing rate. Our main contribution to this constant size setting is that our assumptions allow for hard killing events (a feature observed in several models, such as diffusion processes killed at the boundary of a domain), that we consider general Markov processes (allowing implicit dependence with respect to time and with respect to the past of the process), and do not require continuity of the branching/killing rates with respect to the empirical measure of the process. As is commonly observed in the literature, uniform convergence with respect to time can be obtained assuming additional mixing conditions on the semigroup Q (see the above references and also [START_REF] Del Moral | Exponential mixing properties for time inhomogeneous diffusion processes with killing[END_REF][START_REF] Oçafrain | Convergence of a non-failable mean-field particle system[END_REF][START_REF] Tough | The Fleming-Viot process with McKean-Vlasov dynamics[END_REF] for models with hard killing). We illustrate this property in Section 3.1.

Remark 11. It is also possible to treat branching processes with non-local branching in this setting. Indeed, consider a branching process where, between branching events, particles move according to a Markov process (ξ t ) t≥0 from their point of creation, and when at position x ∈ E, at rate β(x), the particle is removed from the system and replaced by two particles at (random) positions x 1 , x 2 ∈ E, which may be different to x. Particles also undergo (soft) killing at rate κ(x). Denoting by (ψ t ) t≥0 the linear semigroup of this process, in this case, the many-to-one formula is given by

ψ t [g](x) = E x e t 0 b(Xs)-κ(Xs)ds g(X t )1 t<τ ∂ ,
where (X t ) t≥0 behaves like (ξ t ) t≥0 , except that at rate 2b(X t ), the process X t jumps to a new location given by (δ x 1 +δ x 2 )/2. We refer the reader to section 3.2 for an example of a branching process with non-local branching.

Remark 12. The main property of the process ( St , (X i t ) t∈ St ) t≥0 that is used in the proof of Theorem 6 is the identity

E t x Q T -(t+s) f (X 1 t+s )1 t+s<τ 1 = E t Q T -(t+s) f (Y 1 t+s )e -t+s t hu du 1 t+s<τ ∂ 1 t<τ 1 , (8)
where E t is the expectation on a probability space Ω such that (Y i u ) u≥t , i ∈ S0 , are independent copies of X, starting from X i t at time t, and where

h u := j∈ S0 b j ((Y i u ) i∈ S0 ) + κ j ((Y i u ) i∈ S0 ).
This identity is derived from the Markov property, and the assumption on the dynamics of the process ( St , (X i t ) t∈ St ) t≥0 . As such, it is possible to see that Theorem 6 still holds under milder conditions on the construction. Specifically, suppose that there exist a filtered probability space (Ω, F, (F t ) t≥0 , P) and a progressively measurable process ( St , (X i t ) t∈ St ) t≥0 defined on this space, such that the rates b i , κ i , p i , q i are no longer functions, but merely progressively measurable processes, satisfying the equivalents 1 of (2), (3):

b i t -κ i t = b(X i t ) -κ(X i t ) p i t = 0 whenever | St | = 1.
Then the proof of Theorem 6 still holds, provided that condition (8) can be verified in the more general case, and where we interpret E t

x as the expectation conditional on the full sigmaalgebra F t . Of course, such a property is something that may, in certain settings, be obtained through construction.

An example where such a construction might be natural is in a distributed computing environment. In this case, one may look to break the particles into subsets which are handled by separate processors with minimal communication. In this case, one may choose to impose N min -N max -like criterion on each separate processor. Roughly, each processor might be allocated an initial set of particles, which evolve without communication (by only having killing/branching events), so long as the number of particles remains within [N min , N max ]. Since global properties of the particle system are only needed when resampling or selection events occur, the need for communication events in the system, which cause computational bottlenecks, could be minimised by choosing to allow these events provided the number of particles on a given processor does not move outside a specified range. In such a case, the number of particles on a given processor would depend on the allocation of particles to given processors (and which may itself be randomised independently of the history of the particle system). Then the processes b i and κ i would not be simply functions of the particle positions (since we also need to know which particles are on the same processor), but would be expected to be (F t ) t≥0 -progressively measurable, for an appropriate filtration, which may then be larger than the filtration generated by the particle histories.

Applications

Population size dynamics with constrained population size

In this section, we introduce a branching model for a population where the size is constrained between N min ≥ 2 and N max ≥ N min and demonstrate that the results from the previous section may be applied. In this model, we associate with each individual a health state x ∈ E := Z d + , d ≥ 1, which evolves according to the dynamics of a multi-dimensional birth and death process during the life of an individual. The parameters of the pure jump Markovian dynamics are denoted by

β i (x) = transition rate from x to x + e i , δ i (x) = transition rate from x to x -e i ,
1 Technically, we need that the integrals of the respective processes are equal up to modification, that is,

if we define I i,1 t := t 0 (b i s -κ i s ) ds and I i,2 t := t 0 (b i (X i s ) -κ(X i s )
) ds, then we require the processes I i,1 and I i,2 to agree up to a modification. Note that as a consequence of the difference b(•) -κ(•) being bounded, this implies that the processes are also indistinguishable.

where e i ∈ Z d + denotes the vector with 1 in the i th position and 0 everywhere else, β i , δ i : Z d + → [0, +∞) for all i ∈ {1, . . . , d}, and δ i (x) = 0 for all x = (x 1 , . . . , x d ) such that x i = 0 (so that the process cannot leave the set Z d + ). When its health state is x ∈ Z d + , at rate κ(x) an individual dies, where κ : Z d + → [0, +∞), and at rate b(x), it gives birth to an identical individual with equal health state x, where b : Z d + → [0, +∞) is bounded. Note that, in our model, a higher health state may thus be associated with a worse health condition. This corresponds, for instance, to the situation where x is the size of a (harmful) parasite population or the concentration of a virus.

In addition, an external mechanism controls the total size of the individuals population: the total number of individuals is bounded above by N max by removing an individual chosen randomly and uniformly in the population when the population size reaches N max + 1, and the total number of individuals is lower bounded by N min by cloning an individual chosen randomly and uniformly when the population size reaches N min -1 (in the situation where N min = 0, this last mechanism does not operate).

Formally, the state space of our population model is

F := Nmax N =N min [E] N , where [E] N denotes
the unordered N -tuples of E. We represent the elements of F by x = [x 1 , . . . , x N ] and note that repetitions are allowed. The extended infinitesimal generator of the process acting on bounded mesurable functions f :

F → R, is given, for any x = [x 1 , . . . , x N ] ∈ F , by Lf (x) = N k=1 d i=1 β i (x k,i )∆ x,x+e N k,i f + δ i (x k,i )∆ x,x-e N k,i f + N k=1 b(x k )∆ x,x∪[x k ] f + N k=1 κ(x k )∆ x,x\[x k ] f (x) + 1 N =N min N k=1 κ(x k ) 1 N -1 N =1, =k ∆ x\[x k ],x∪[x ]\[x k ] f + 1 N =Nmax N k=1 b(x k ) 1 N + 1 N =1 (1 + 1 =k )∆ x∪[x k ],x∪[x k ]\[x ] f, where • x k,i is the i th component of x k , • x + e N k,i is x but with x k,i replaced by x k,i + 1, • x -e N k,i is x but with x k,i replaced by x k,i -1, • for all x, y, ∆ x,y f := f (y)1 y∈F -f (x)1 x∈F , • for all y ∈ E, x ∪ [y] = [x 1 , . . . , x N , y], • x \ [x k ] = [x 1 , . . . , x k-1 , x k+1 , . . . , x N ].
In a similar manner to the previous section, we denote by A T the number of times an individual is duplicated by the size constraining mechanism before time T , and by B T the number of times an individual is removed from the population by the size constraining mechanism. We also set (X(t)) t≥0 to be the population process, (|X(t)|) t≥0 its size, and (X t ) t∈≥0 a multidimensional birth and death process with (state dependent) transition rates β and δ. The notation P x (respectively P x ) is used for the law of the population process X starting from x ∈ F (respectively for the law of X starting from x ∈ E). We denote as usual by E x and E x the associated expectations.

The first difficulty is to ensure that the constrained branching model defined above does not explode in finite time, where explosion is interpreted in the sense that infinitely many jumps of the system happen in finite time. Indeed, since κ is not assumed to be bounded, the rate at which new individuals are added to the system is not bounded. The following result provides a sufficient criterion for the non-explosion of the process but before stating it, let us first introduce some definitions that will be used in what follows.

We say that two functions V and κ are co-monotone if (κ(x) -κ(y))(V (x) -V (y)) ≥ 0 for all x, y ∈ E. We say that V and κ are almost co-monotone if there exists a function κ co-monotone with V such that κ -κ is bounded. We say that V goes to infinity at infinity if the set of x ∈ E such that V (x) < c is finite for all c > 0.

Proposition 13. Assume that there exists a function V : E → [1, +∞) almost co-monotone with κ, such that V tends to infinity at infinity, and such that, for all x ∈ E,

d i=1 [β i (x)(V (x + e i ) -V (x)) + δ i (x)(V (x -e i ) -V (x))] ≤ CV (x), (9)
for some constant C ∈ R. Then the process with transition rates given by L is non-explosive.

Remark 14. The proof of this non-explosion result makes use of the extended infinitesimal generator of the particle system and of Dynkin's formula. As a consequence, its proof can be extended without too much effort to more general continuous time Markov processes, by replacing the left hand term in (9) by the extended infinitesimal generator of the more general process. One difficulty can appear though when the process is not a simple pure jump process on a discrete state space: to prove in a general context that L (adapted to fit the general individual dynamics of the particles) is indeed the extended infinitesimal generator of the particle system can be challenging.

Remark 15. More complex resampling/selection mechanisms could be considered with only slight modifications of the proof of Proposition 13. For instance, the main line of the proof applies as well if N max = +∞, or if we assume that, immediately after each branching time and with probability 1 + 1

N max + 1 - 1 N
, a uniformly chosen individual is removed from the population, and/or that, immediately after each killing time and with probability

N min -1 N , a
uniformly chosen individual is duplicated (where N is the population size just after the event).

Note that these choices of parameters may be used to model an increasing competition between individuals when the size of the population increases.

Proof of Proposition 13. We define

V :F → [1, +∞) :x = [x 1 , . . . , x N ] → N k=1 V (x k ),
so that V goes to infinity at infinity and for

x = [x 1 , . . . , x N ] ∈ F , L V (x) = N k=1 d i=1 [β i (x k )(V (x k + e i ) -V (x k )) + δ i (x k )(V (x k -e i ) -V (x k ))] + N k=1 (b(x k ) -κ(x k )) V (x k ) + 1 N =N min N k=1 κ(x k ) 1 N -1 N =1, =k V (x l ) -1 N =Nmax N k=1 b(x k ) 1 N + 1 N =1 (1 + 1 =k )V (x ).
Now let κ be co-monotone with V such that κ = κ -κ is bounded. Then, using the fact that b is bounded, we deduce that

L V (x) ≤ (C + b ∞ )V (x) + 1 N =N min    1 N -1 N k=1 κ(x k ) N =1 =k V (x ) - N k=1 κ(x k )V (x k )    = (C + b ∞ )V (x) + 1 N =N min N -1 N k=1 κ(x k ) N =1 V (x ) -N N k=1 κ(x k )V (x k ) = (C + b ∞ )V (x) + N 2 1 N =N min N -1 N k=1 κ(x k ) N N k=1 V (x k ) N - N k=1 κ(x k )V (x k ) N ≤ (C + b ∞ + 4 κ ∞ )V (x) + N 2 1 N =N min N -1 N k=1 κ (x k ) N N k=1 V (x k ) N - N k=1 κ (x k )V (x k ) N ≤ (C + b ∞ + 4 κ ∞ )V (x), (10)
where we used the fact that κ and V are co-monotone and the FKG inequality for the last inequality. For all n ≥ 1, let

U n = {x ∈ F, max k∈{1,...,|x|}
x k > n} and define the random variable

θ n := inf{t ≥ 0, (X 1 t , . . . , X N t ) ∈ U n }.
Note that θ n is a stopping time since the constrained branching process is càdlàg. Moreover, denoting by C n the set of points in F that can be reached in one jump from F \ U n , we observe that C n is bounded and hence that V (X t∧θn ) is almost surely uniformly bounded in t ≥ 0. In particular, we deduce, using Dynkin's formula and the inequality [START_REF] Champagnat | Exponential convergence to quasi-stationary distribution and Q-process[END_REF], that, for all n ≥ 0, all t ≥ 0 and all x ∈ F ,

E x ( V (X t∧θn )) ≤ V (x) + (C + b ∞ + 4 κ ∞ ) t 0 E x ( V (X s∧θn )) ds.
Hence, using Grönwall's inequality, we deduce that

E x ( V (X t∧θn )) ≤ V (x)e t(C+ b ∞+4 κ ∞) . Since {X t∧θn ∈ U n , θ n ≤ t} = {θ n ≤ t} up to a negligible event, we conclude that inf y∈Un V (y)P x (θ n ≤ t) ≤ V (x) e t(C+ b ∞+4 κ ∞) .
Denote by τ ∞ the explosion time of the process. Since the jump rate of the process is bounded up to time θ n , we deduce that θ n ≤ τ ∞ almost surely. Hence

P x (τ ∞ ≤ t) ≤ V (x)e t(C+ b ∞+4 κ ∞) inf y∈Un V (y) ----→ n→+∞ 0,
so that τ ∞ = +∞ almost surely. This concludes the proof.

Since the process defined in this section is a particular instance of the process described in Section 2.1, the following corollary is a direct consequence of Proposition 13 and Theorem 6.

Corollary 16. Under the assumptions of Proposition 13, we have, for all T > 0, all bounded measurable functions f :

Z d + → R and all initial position X(0) ∈ F , E X(0)   Π A T Π B T |X(t)| k=1 f (X(t) k )   = |X(0)| k=1 E X(0) k f (X t ) exp t 0 (b(X s ) -κ(X s )) ds
where

Π A T = N min -1 N min A T and Π B T = N max + 1 N max B T .
In addition,

|X(0)| k=1 E X(0) k f (X t ) exp t 0 (b(X s ) -κ(X s )) ds |X(0)| k=1 E X(0) k exp t 0 (b(X s ) -κ(X s )) ds - |X(t)| k=1 f (X(t) k ) |X(t)| 2 ≤ C T f ∞ 1 |X(0)| |X(0)| k=1 E X(0) k exp t 0 (b(X s ) -κ(X s )) ds 1 |X(0)| .
Our aim is now to study the long-time behaviour of the (normalised) empirical distribution of the health parameter across the population. We assume that there exist constants

C 1 > 0, C 2 > 0, η 0 > 0, n 0 ≥ 0 such that, for all x ∈ Z d + satisfying |x| ≥ n 0 , d i=1 (δ i (x) -β i (x)) ≥ C 1 |x| 1+η 0 (11) d i=1 (δ i (x) -β i (x)) ≥ C 2 ( b ∞ -b(x) + κ(x))|x|. ( 12
)
Proposition 17. Assume that (11) and (12) hold true and that X is irreducible. Then there exists a probability distribution ν X on Z d + such that, for all probability measures µ on Z d + and all bounded measurable functions f on Z d + ,

E µ f (X t ) exp t 0 b(X s ) -κ(X s ) ds E µ exp t 0 b(X s ) -κ(X s ) ds -ν X (f ) ≤ C X e -γ X t f ∞ , (13)
for some constants C X ≥ 1 and γ X > 0. In addition, there exists a bounded function η X :

Z d + → (0, +∞) and λ X > -b ∞ such that for all probability measures µ on Z d + , E µ η X (X t ) exp t 0 b(X s ) -κ(X s ) ds = e -λ X t µ(η X ), ∀t ≥ 0.
Proof of Proposition 17. First observe that, for all bounded measurable functions f :

Z d + → R, E µ f (X t ) exp t 0 b(X s ) -κ(X s ) ds = e b ∞t E b µ f (X b t )1 t<τ b ∂ ,
where P b µ and E b µ denote the law and associated expectation of the birth and death process X b with birth and death parameters β and δ, killed at rate κ b := b ∞ -b + κ at a time denoted by τ b ∂ . We now prove that there exist ν X , C X > 0 and γ X > 0 such that, for all bounded measurable functions f :

Z d + → R E b µ f (X b t ) t < τ b ∂ -ν X (f ) ≤ C x e -γ X t f ∞ , (14)
which is equivalent to [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF].

In order to prove ( 14), we check that the assumptions of [11, Theorem 2.4] hold true. Since the proof is similar to Theorem 4.1 in the above reference, we only recall the main ingredients. Set W (x) = |x| for all x ∈ Z d + . Assumption [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF] implies that, for all x ∈ Z p + large enough,

d i=1 [β i (x)(W (x + e i ) -W (x)) + δ i (x)(W (x -e i ) -W (x))] ≤ -C 1 W (x) 1+η 0 .
This classically entails that the birth and death process on Z d + with birth rates β and death rates δ (without absorption) is non-explosive and comes down from infinity. This implies that the process is regularly absorbed (as defined in (2.1), p.55 of [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF]) and that Assumption 3 in [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF] is satisfied.

Moreover, choosing appropriately a 1 > 1, a 2 > 1 and ε > 0, and setting

V (x) = |x| k=1 1 k a 1 and ϕ(x) = +∞ k=|x|+1 1 k a 2 ,
one deduces, using the same computations as in Section 4 in [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF], that Assumption 1 of [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF] is satisfied. The only difference is that the infinitesimal generator that one needs to consider in the present case is

Lf (x) = d i=1 [β i (x)(f (x + e i ) -f (x)) + δ i (x)(f (x -e i ) -f (x))] -κ b (x)f (x).
Finally, Assumption 2 in the above reference is immediate since the process is irreducible in a locally finite state space. This shows that [START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed markov processes[END_REF]Theorem 2.4] applies to X b , which concludes the proof of [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF].

The last assertion of the proposition follows from [10, Proposition 2.3]: it implies that there exists λ 0 > 0 and a bounded positive function η X such that

E b µ (η X (X b t )) = e -λ 0 t η X .
Multiplying by e b ∞t and setting λ X = λ 0 -b ∞ concludes the proof.

We consider now the situation where κ is bounded and where N min is large, and prove that the normalised empirical distribution

m T m T (1 E )
of the population health's parameter converges toward ν X when N min , N max and T go to infinity.

Proposition 18. Assume that the conditions of Proposition 17 hold and that κ is bounded. Then, denoting by m 0 the initial empirical measure of the population's health parameter and setting m 0 = m 0 /m 0 (1 E ), we have Proposition 17, and C, c are from Theorem 6. In particular,

sup T ∈[0,+∞) m T (f ) m T (1 E ) - E m 0 f (X t ) exp t 0 (b(X s ) -κ(X s )) ds E m 0 exp t 0 (b(X s ) -κ(X s )) ds 2 ≤ (2C X +C)N -α/2 min f ∞ , where α = γ X /((c + 1) b ∞ + κ ∞ + γ X ) ∈ (0, 1), C X , γ X are from
m T (f ) m T (1 E ) -ν X (f ) 2 ≤ (2C X + C)N -α/2 min f ∞ + C X e -γ X T f ∞ .
Remark 19. This result can be adapted to include the case where X is a Markov process with hard killing at a boundary when the property obtained in Proposition 17 is known to hold true (see for instance [START_REF] Del Moral | On the stability of nonlinear Feynman-Kac semigroups[END_REF][START_REF] Champagnat | Exponential convergence to quasi-stationary distribution and Q-process[END_REF][START_REF] Champagnat | Criteria for exponential convergence to quasi-stationary distributions and applications to multi-dimensional diffusions[END_REF][START_REF] Horton | Stochastic methods for the neutron transport equation I: Linear semigroup asymptotics[END_REF][START_REF] Del Moral | Exponential mixing properties for time inhomogeneous diffusion processes with killing[END_REF]), provided one can prove that the particle system does not degenerate toward the boundary. However this can be challenging, see for instance [START_REF] Grigorescu | Ergodic properties of multidimensional Brownian motion with rebirth[END_REF][START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF][START_REF] Bieniek | Extinction of Fleming-Viot-type particle systems with strong drift[END_REF].

Proof of Proposition 18. We assume without loss of generality that f ∞ = 1 and that N min > N max (otherwise the result is trivial). We define the semigroup Q by

Q t f (x) = E x f (X t ) exp t 0 (b(X s ) -κ(X s )) ds .
We have, by Proposition 17 applied first to

µ = m 0 Q T -t /m 0 Q T -t 1 E and then to m T -t = m T -t /m T -t 1 E ,
and by Theorem 6, for all 0 ≤ t ≤ T ,

m 0 Q T f m 0 Q T 1 E - m T (f ) m T (1 E ) 2 ≤ m 0 Q T -t Q t f m 0 Q T -t Q t 1 E - m T -t (Q t f ) m T -t (Q t 1 E ) 2 + m T -t (Q t f ) m T -t (Q t 1 E ) - m T (f ) m T (1 E ) 2 ≤ 2C X e -γ X t + C exp(c b ∞ t) Q t 1 E ∞ m T -t Q t 1 E 1 √ m T -t 1 E 2 , (15)
where

Q t 1 E ∞ m T -t Q t 1 E 1 √ m T -t 1 E 2 ≤ exp( b ∞ t) exp(-κ ∞ t) 1 √ N min . Set t 0 = ln √ N min (c + 1) b ∞ + κ ∞ + γ X ,
so that e -γ X t 0 = e ((c+1) b ∞+ κ ∞)t0 / N min . If T ≥ t 0 , then choosing t = t 0 in (15), we get

m 0 Q T f m 0 Q T 1 E - m T (f ) m T (1 E ) 2 ≤ (2C X + C)e -γ X t 0 = (2C X + C)N -α/2 min , with α = γ X (c + 1) b ∞ + κ ∞ + γ X
. For T ≤ t 0 , Theorem 6 implies that

m 0 Q T f m 0 Q T 1 E - m T (f ) m T (1 E ) 2 ≤ C exp ((c + 1) b ∞ t 0 + κ ∞ t 0 ) 1 √ N min = CN -α/2 min
which concludes the proof of the first inequality in Proposition 18. The second inequality is then a straightforward consequence of the first inequality and of Proposition 17.

Piecewise deterministic Markov processes

Due to Assumption 1 in Section 2, the results presented thus far do not hold for the case where (X t ) t≥0 is a piecewise deterministic Markov process (PDMP) and E is a bounded domain with absorbing boundary conditions, since two independent copies of the process may hit ∂E at the same time. In this section, we consider a particular example, namely the neutron transport equation (NTE), in order to discuss a way to avoid this problem.

The NTE is a balance equation that describes the behaviour of neutrons in a fissile medium such as a nuclear reactor. In such systems neutrons move in straight lines with a fixed speed until they either come into contact with the boundary of the reactor, at which point they are absorbed, or they collide with the nucleus of an atom. When the latter occurs, either the neutron undergoes a scattering event where the particle bounces off the nucleus and continues its motion but with a new velocity, or a fission event occurs where the collision causes new neutrons to be produced with identical spatial positions, but potentially different velocities.

In this setting, we have E = D × V where D ⊂ R 3 is the spatial domain, which we assume to be open and bounded with smooth boundary, ∂D, and V = {υ ∈ R 3 : v min < |υ| < v max } is the velocity domain, where 0 < v min ≤ v max < ∞. The NTE can then be stated as follows:

∂ ∂t ψ t (r, υ) = υ • ∇ψ t (r, υ) -(σ s (r, υ) + σ f (r, υ))ψ t (r, υ) + σ s (r, υ) V ψ t (r, υ )π s (r, υ, υ )dυ + σ f (r, υ) V ψ t (r, υ )π f (r, υ, υ )dυ , (16)
where σ s (r, υ) : the rate at which scattering occurs from incoming velocity υ at position r, σ f (r, υ) : the rate at which fission occurs from incoming velocity υ at position r, π s (r, υ, υ ) : probability density that an incoming velocity υ at position r scatters to an outgoing velocity, with probability υ satisfying V π s (r, υ, υ )dυ = 1, and π f (r, υ, υ ) : density of expected neutron yield at velocity υ from fission with incoming velocity υ satisfying V π f (r, υ, υ )dυ < ∞.

Moreover, we impose the following initial and boundary conditions

ψ 0 (r, υ) = g(r, υ) for r ∈ D, υ ∈ V, ψ t (r, υ) = 0 for r ∈ ∂D if υ • n r > 0,    ( 17
)
where n r is the outward facing normal of D at r ∈ ∂D and g :

D × V → [0, ∞) is a bounded, measurable function.
Under the assumptions, (A1) σ s , σ f , π s and π f are uniformly bounded away from infinity.

(A2) inf r∈D,υ,υ ∈V (σ s (r, υ)π s (r, υ, υ ) + σ f (r, υ)π f (r, υ, υ )) > 0,
it was shown in [START_REF] Cox | Multi-species neutron transport equation[END_REF][START_REF] Horton | Stochastic methods for the neutron transport equation I: Linear semigroup asymptotics[END_REF] that the NTE can be modelled using a weighted PDMP, also known as the neutron random walk (NRW). More precisely, setting

α(r, υ) = σ s (r, υ) + σ f (r, υ) V π f (r, υ, υ )dυ (18) π(r, υ, υ ) = (α(r, υ)) -1 [σ s (r, υ)π s (r, υ, υ ) + σ f (r, υ)π f (r, υ, υ )] , (19)
β(r, υ) = α(r, υ) -σ s (r, υ) -σ f (r, υ), (20)
we define the NRW ((R t , Υ t ) t≥0 , P (r,υ) ) as follows. From an initial configuration (r, υ) ∈ D ×V , the particle will propagate linearly until it is either absorbed at the boundary of the domain or, at rate α the process scatters off a nucleus and a new velocity is chosen according to π.

Then, for a bounded, measurable function g, and (r, υ)

∈ D × V ψ t [g](r, υ) := E (r,υ) e t 0 β(Rs,Υs)ds g(R t , Υ t )1 (t<τ D ) ,
solves ( 16)- [START_REF] Dawson | Measure-valued Markov processes[END_REF], where τ D is the first exit time of the NRW from D.

It was also shown in [START_REF] Horton | Stochastic methods for the neutron transport equation I: Linear semigroup asymptotics[END_REF] that if (A1) and (A2) are satisfied, and that if there exists an ε > 0 such that (B1) D ε := {r ∈ D : inf y∈∂D |r -y| > εv max } is non-empty and connected, (B2) there exist 0 < s ε < t ε and γ > 0 such that, for all r ∈ D\D ε , there exists K r ⊂ V measurable such that Vol(K r ) ≥ γ > 0 and for all υ ∈ K r , r + υs ∈ D ε for every s ∈ [s ε , t ε ] and r + υs / ∈ ∂D for all s ∈ [0, s ε ],

the semigroup (ψ t ) t≥0 exhibits the following Perron Frobenius behaviour.

Theorem ( [START_REF] Horton | Stochastic methods for the neutron transport equation I: Linear semigroup asymptotics[END_REF]). There exists a λ * ∈ R, a positive right eigenfunction ϕ ∈ L + ∞ (D × V ) and a left eigenmeasure which is absolutely continuous with respect to Lebesgue measure on D × V with density φ ∈ L + ∞ (D × V ), both having associated eigenvalue e λ * t , and such that ϕ (resp. φ) is uniformly (resp. a.e. uniformly) bounded away from zero on each compactly embedded subset of D × V . In particular, for all

g ∈ L + ∞ (D × V ), φ, ψ t [g] = e λ * t φ, g (resp. ψ t [ϕ] = e λ * t ϕ), t ≥ 0. ( 21
)
Moreover, there exists ε > 0 such that

sup g∈L + ∞ (D×V ): g ∞≤1 e -λ * t ϕ -1 ψ t [g] -φ, g ∞ = O(e -εt ), t ≥ 0. ( 22
)
This theorem therefore provides a way to build Monte Carlo simulations to estimate the eigenvalue λ * and the eigenfunctions ϕ and φ via the NRW. We refer the reader to [START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF] for further details. However, although simulating a single weighted path has advantages over simulating an entire tree of neutrons, the transience of this process means that many of the particles exit the domain relatively quickly and therefore, a large number of simulations are required in order to obtain information about the system. This also leads to problems with the variance of the estimators. In order to deal with this problem, the notion of 'h-transform' was developed in [START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF], where one biases the NRW to prevent it from exiting the domain.

More precisely, let h be a bounded positive function on D × V such that h = 0 on ∂D and define the following change of measure,

dP h (r,υ) dP (r,υ) σ((Rs,Υs),s≤t) := exp - t 0 Jh(R s , Υ s ) h(R s , Υ s ) ds Nt i=1 h(R T i , Υ T i ) h(R T i , Υ T i-1 ) 1 (t<τ D ) , (23)
where

Jg(r, υ) = α(r, υ) V [g(r, υ ) -g(r, υ)]π(r, υ, υ )dυ , (24)
(T i ) i≥0 are the scatter times of the NRW with T 0 = 0, and N t = sup{i : T i ≤ t}.

Then, ((R, Υ), P h ) also defines a NRW but with scattering operator

J h g(r, υ) = α(r, υ) V [g(r, υ ) -g(r, υ)] h(r, υ ) h(r, υ) π(r, υ, υ )dυ . ( 25
)
The idea is that the factor h(r, υ )/h(r, υ) in the above integral forces the NRW to scatter when it is approaching the boundary, thus preventing it from being killed. Thus, if we can show that this process satisfies Assumption 1, we may use the BBMMI process in conjunction with the h-transformed process to estimate the leading eigentriple (λ * , ϕ, φ), for example. Indeed, it is straightforward to show that the leading eigentriple of the h-transformed process is given by (λ * , ϕ/h, h φ) and so we may use the BBMMI process to estimate this triple, and hence obtain an estimate for the original quantities (λ * , ϕ, φ). The rest of this section is dedicated to verifying Assumption 1 for the transformed process. We refer the reader to Section 3.3 for the numerical aspects.

To this end, fix δ = 1/(2 α ∞ ) and let φ : [0, +∞) → [0, δ] be a regular non-decreasing function such that φ(x) = x for all x ∈ [0, δ/2] and φ(x) = δ for all x ≥ δ. We set h(r, υ) = φ (κ r,υ ) .

Using the following result, one can express expectations with respect to P as Feynman-Kac formulas related to P h . In what follows, we say that (r, υ) is in a vicinity of ∂D if h(r, υ) < δ/2. Lemma 20. Suppose (A1), (A2) hold. A process with law P h (r,v) does not hit the boundary with probability one. Moreover, setting b = Lh h + and κ = Lh h - , the function b is bounded on D × V and κ is locally bounded. Let P κ r,υ be the law of a process with law P h r,υ with additional soft killing at rate κ, then we have

E (r,υ) f (R t , Υ t ) 1 (t<τ D ) = h(r, υ) E κ (r,υ) exp t 0 b(R s , Υ s ) ds f (R t , Υ t ) h(R t , Υ t ) 1 (t<τκ) , (26)
where τ κ denotes the soft killing time of the process with law P κ r,υ . Proof. We have

ψ t [g](r, υ) = E h (r,υ) exp t 0 Jh(R s , Υ s ) h(R s , Υ s ) + β(R s , Υ s )ds Nt i=1 h(R T i , Υ T i-1 ) h(R T i , Υ T i ) g(R t , Υ t )1 (t<τ D ) (27)
Note that inf r∈D,υ∈V

lim s→κ D r,υ |υ|(κ D r,υ -s) h(r + υs, υ) > 0, (28)
where κ D r,υ := inf{t > 0 : r + υt ∈ D}. Hence, Theorem 7.1 in [START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF] implies that, under P h (r,υ) , r ∈ D, υ ∈ V , (R, Υ) does not hit the boundary with probability one.

Recall that the infinitesimal generator associated to ((R t , Υ t ) t≥0 , P) is given by

Lf (r, υ) = Tf + α(r, υ) V (f (r, υ ) -f (r, υ)) π(r, υ, υ ), dυ,
where T is defined for regular functions by Tf (r, υ) = υ • ∇ r f (r, υ) and, more generally and when the limit exists, by

Tf (r, υ) = lim ε→0 + f (r + ευ) -f (r, υ) ε .
In particular, the conditions of Remark 7.1 of [START_REF] Cox | Monte-Carlo methods for the neutron transport equation[END_REF] are satisfied and hence dP h (r,υ)

dP (r,υ) σ((Rs,Υs),s≤t) = exp - t 0 Lh(R s , Υ s ) h(R s , Υ s ) ds h(R t , Υ t ) h(r, υ) 1 (t<τ D ) . (29)
Finally, we obtain, for all measurable functions f :

D × V → [0, +∞), E (r,υ) f (R t , Υ t ) 1 (t<τ D ) = h(r, υ) E h (r,υ) exp t 0 Lh(R s , Υ s ) h(R s , Υ s ) ds f (R t , Υ t ) h(R t , Υ t ) ( 30
)
where have used the fact that 1 (t<τ D ) = 1, P h (r,υ) -almost surely, deduced from [START_REF] Etheridge | Some Mathematical Models from Population Genetics: École D' Été de[END_REF]. Now, introducing b = Lh h + and κ = Lh h - , and considering the law P κ r,υ of a process with law P h r,υ with soft killing at rate κ, we deduce that

E (r,υ) f (R t , Υ t ) 1 (t<τ D ) = h(r, υ) E κ (r,υ) exp t 0 b(R s , Υ s ) ds f (R t , Υ t ) h(R t , Υ t ) 1 (t<τκ) , (31)
where τ κ denotes the soft killing time of the process (with law P h r,υ with soft killing at rate κ). Now, we observe that, in a vicinity of ∂D,

Lh(r, υ) = -1 + α(r, υ) V (h(r, υ ) -h(r, υ)) π(r, υ, υ ) dυ (32) ≤ -1 + α ∞ δ = -1/2
and hence that sup

(r,υ)∈D×V Lh(r, υ) h(r, υ) < +∞.
This implies that b is bounded. The fact that κ is locally bounded is an immediate consequence of the regularity of h.

Our aim is now to apply the results of Section 2.2 to the Feynman-Kac expression ( 26) in order to obtain information about P. Assumption 1 is clearly satisfied by a process with law P κ (r,υ) , and it only remains to check that Assumption 2 holds true in order to apply Theorem 6. This is the purpose of the following result. In what follows, we refer to a bounded BBMMI as any BBMMI whose selection mechanism imposed by the parameters p i results in a size constrained process, i.e. there is a constant N max such that Nt ≤ N max for all t ≥ 0. Proposition 21. Under the assumptions (A1), (A2), (B1) and (B2), any bounded BBMMI driven by P κ (r,υ) satisfies Assumption 2.

Proof. We denote by σ 1 < σ 2 < ... the sequence of times at which events occur for a given particle in the system (each time is either a scattering, a branching or a soft killing). Note that each time σ n is well defined since, almost surely, the number of event occurring before a time T goes to infinity when T → +∞. Since one can use the strong Markov property at the birth time of the given particle, we assume without loss of generality that the particle is already alive at time 0. From the expression of the scattering operator J h in [START_REF] Del Moral | Exponential mixing properties for time inhomogeneous diffusion processes with killing[END_REF], one observes that the scattering rate toward a direction in a set V ⊂ V is given by V α(r, υ) h(r, υ ) h(r, υ) π(r, υ, υ )dυ . But, under the assumptions (B1-2), there exists ε > 0, β > 0 and η > 0 such that, for any point r ∈ D, there exists V r with r

+ [0, β]V r ⊂ D and Leb d ((r + [0, β]V r ) ∩ {r ∈ D, d(r , ∂D) ≥ ε}) > η,
where Leb d denotes the d-dimensional Lebesgue measure. In particular, using the fact that h(r, υ ) is lower bounded on V r (uniformly over r ∈ D). We deduce that there exists a constant s > 0 such that the scattering rate toward a direction in V r is lower bounded by s/h(r, υ), uniformly in (r, υ) ∈ D × V . Using [START_REF] Grigorescu | Hydrodynamic limit for a Fleming-Viot type system[END_REF] to bound Lh from below, there exists a constant c such that the total rate at which a particle undergoes an event (either a scattering, a branching or a (soft) killing) is upper bounded by c/h(r, υ), uniformly in (r, υ) ∈ D × V . Hence each event has a probability greater than s/c to be a scattering event toward a direction in V r , independently of the past of the process, and in particular independently of the time event σ n : formally, for all n ≥ 0,

P κ σ n is a scattering event and (R σn , Υ σn ) ∈ {R σn-} × V R σn-| σ n , R σn-, Υ σn-≥ s /c > 0,
where (R, Υ) denotes the position and direction of the given particle. Now, using the fact that the killing rate and the scattering rate are uniformly bounded in {r ∈ D, d(r , ∂D) ≥ ε} and the fact that the total branching rate of the system is uniformly bounded, we deduce that there exists a constant p > 0 and a time T > 0 such that, for all r ∈ D and υ ∈ V r , P(σ 1 ≥ T ) ≥ p.

Hence, we deduce from the strong Markov property that

P κ (σ n+1 < T ) ≤ E κ 1 σn<T P κ Xσ n (σ 1 < T ) ≤ E κ 1 σn<T 1 -1 Υσ n ∈V Rσ n P κ Xσ n (σ 1 < T ) ≤ E κ 1 σn<T 1 -1 Υσ n ∈V Rσ n p ≤ P κ (σ n < T )(1 -s /c p).
This shows that the number of events occurring before time T is stochastically dominated by a geometric random variable with parameter s /c p > 0, so that it is finite almost surely and that its expectation is bounded by a constant ā. Since this is true for all particles, we deduce that total number of events has an expectation bounded by N max ā. Finally, observing that this bound does not depend on the initial position, we deduce using the Markov property at times T, 2T, . . . that the number of event occurring before a time horizon nT is bounded by kN max ā, and hence that it is finite almost surely. This concludes the proof of Proposition 21.

Numerical properties of the constrained branching process

The aim of this section is to discuss some of the numerical properties of the N min -N max process defined in Section 2.1, see Remark 4. In Section 3.3.1, we first make comparisons between this algorithm and the fixed size Moran type algorithm studied in [START_REF] Moral | Feynman-Kac formulae. Probability and its Applications[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF]. We will often use the abbreviation FV IPS to refer to the latter, in reference to the name Fleming Viot interacting particle system, which has been used in several instances in the literature. In Section 3.3.2, we introduce a filtering method to deal with numerical discrepancies occurring in situations where Π A T Π B T has high variance.

Comparison with the FV IPS

We recall that the FV IPS simulates a set of N sub-Markov processes, until the first time one of the particles dies, at which point a resampling event occurs by uniformly selecting one of the remaining N -1 particles and duplicating it. Note that the FV IPS is a particular instance of the N min -N max IPS, with a null branching rate and N min = N max = N . In particular, Theorem 6 applies to the FV IPS (see [START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF] for an analogous result in the FV case).

In some situations, it is possible to use the FV IPS to approximate quantities associated with branching processes. Indeed, it is well known that for a branching process X = (X t ) t≥0 , under fairly weak assumptions, the following many-to-one lemma holds

ψ t [f ](x) := E δx [ f, X t ] = E x e t 0 β(Ys)ds f (Y t ) ,
where β is determined by the branching rate and mean offspring of X, and the law of Y is uniform amongst all paths of X. In the case where β is uniformly bounded above, one may compensate the right-hand side above by its maximum, β, to introduce a penalised process:

ψ k t [f ](x) := e -βt ψ t [f ](x)
. The penalisation has the effect of killing Y at an additional rate β -β, thus yielding a sub-Markov process with semigroup ψ k t . This then allows one to use the FV IPS to obtain estimates of this semigroup and, in some situations, of its associated quasistationary distribution, for example. In general, the use of the FV IPS for the approximation of quantities related to some non-conservative semigroups is widespread, see for instance [START_REF] Burdzy | Configurational transition in a Fleming-Viot-type model and probabilistic interpretation of Laplacian eigenfunctions[END_REF][START_REF] Moral | Feynman-Kac formulae. Probability and its Applications[END_REF][START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF][START_REF] Grigorescu | Hydrodynamic limit for a Fleming-Viot type system[END_REF][START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF][START_REF] Groisman | Simulation of quasi-stationary distributions on countable spaces[END_REF][START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF].

In this section, we show on a simple example that the N min -N max algorithm allows one to overcome some of the limitations of the FV IPS. For simplicity, we only consider the situation where N min = N max .

A simple branching birth and death process. We consider a birth and death process on E = {1, . . . , M } for some parameter M ∈ N. For a particle occupying state x, one of the following things may occur:

• the particle jumps with rate x 2 , to state 1 ∨ (x -1) with probability x x + 1 and to state

(x + 1) ∧ M with probability 1 x + 1 ,
• at rate b(x) = x, a new particle is produced at the same site, which will continue independently the same behaviour as the original particle.

In this situation where the state space is finite and irreducible, it is well known that the semigroup (ψ t ) t≥0 admits a limiting distribution up to an exponential rescaling (see for instance [START_REF] Darroch | On quasi-stationary distributions in absorbing continuoustime finite Markov chains[END_REF]): there exist a constant λ M ∈ R, a positive function η M : {1, . . . , M } → (0, +∞) and a probability measure ν M on {1, . . . , M } such that

sup f :{1,...,M }, f ∞≤1 e -λ M t ψ t [f ](x) -η M (x)ν M (f ) ≤ C M e -γ M t , t ≥ 0,
for some positive constants C M , γ M . The associated birth and death process with killing evolves on {1, . . . , M } with the same jump rates, without branching, and with the additional killing rate κ(x) = M -x. Denoting as above by ψ k the associated sub-Markov semigroup, we have ψ k t = e -M t ψ t and hence sup

f :{1,...,M }, f ∞≤1 e -(λ M -M )t ψ k t [f ](x) -η M (x)ν M (f ) ≤ C M e -γ M t .
Theoretical convergence of the N min -N max and of the FV IPSs. In what follows, m M T refers to the empirical measure of the N min -N max particle system at time T , associated to the above branching birth and death process, and µ M T refers to the empirical measure of the FV IPS with N := N min = N max particles at time T associated to the above killed birth and death process.

In particular, Proposition 18 can be adapted to our setting and hence, in the situation where

N min = N max = N , m M T (f ) N -ν M (f ) 2 ≤ (2C X + C) f ∞ N α/2 + C X e -γ X T f ∞ .
for some α > 0. Using the fact that the N min -N max process is ergodic, we deduce that its normalised empirical stationary distribution X M N (the normalised empirical distribution of a random vector of particles distributed according to the stationary distribution of the N min -N max process) satisfies 1: For each value of M ∈ {10, 100, 1000, +∞}, we display the bias of the estimator θN , its standard deviation, and the number of events, for the N min -N max algorithm (for which θN = X M N (f M )) and for the FV IPS (for which θN = Y M N (f M )) with N = N min = N max = 10 particles. Note that the FV IPS is not defined when M = +∞. for all function f : {1, . . . , M } → ν M . Similarly, the normalised empirical stationary distribution Y M N of the FV particle system satisfies

X M N (f ) L 2 ----→ N →+∞ ν M (f ), N = 10 |E( θN ) -ν M (f M )| Std( θN ) (A T + B T )/T M =
Y M N (f ) L 2 ----→ N →+∞ ν M (f ).
These theoretical convergence results ensure that both algorithms can be used to approximate the limiting distribution ν M . It remains to compare their performances via numerical simulations.

Numerical comparison of the N min -N max and the FV IPSs. We now consider the problem of numerically approximating the mean of ν M , that is ν

M (f M ) with f M (x) = x for all x ∈ {1, . . . , M }. The estimator is θN = X M N (f M ) (resp. θN = Y M N (f M ))
for the N min -N max (resp. FV) IPS, and we compare three important metrics determining the performance of the algorithms:

• the bias of the estimator, |E( θN ) -ν M (f M )|,
• the standard deviation of the estimator, Std( θN ),

• the number of interaction events, (A T + B T )/T , which corresponds to the number B T of selection events in the N min -N max algorithm and the number A T of resampling events in FV IPS, per unit of time.

The last metric (A T + B T )/T is particularly important since it entails that the total number of interaction events grows (at least) linearly in time.

For N = 10 and then N = 100, we let the upper bound of the state space M vary in {10, 100, 1000, +∞}, and present the results in Table 1 for N = 10, and in Table 2 for N = 100. For each metric, a lower number indicates that the algorithm performs better. For each value of M ∈ {10, 100, 1000, +∞}, we display the bias of the estimator θN , its standard deviation, and the number of events, for the N min -N max algorithm (for which θN = X M N (f M )) and for the FV IPS (for which θN = Y M N (f M )) with N = N min = N max = 100 particles. Note that the FV IPS is not defined when M = +∞.

N = 100 |E( θN ) -ν M (f M )| Std( θN ) (A T + B T )/T M =
We can see that, for this particular example, the N min -N max algorithm performs better in all situations, with respect to all metrics. Moreover, the bias of the FV approximation increases when M increases, while the variance of the estimator also increases and the computation cost increases. Without surprise, increasing M does not change the dynamic of the N min -N max algorithm, since very few particles, if any, reach the boundary M even when M = 10. Further, this suggests that it should be possible to extend the N min -N max algorithm to the case of unbounded branching rates, and this is supported by the numerical simulations when M = +∞ (for which the FV IPS is of course not defined), although it is not covered by our main result. We intend to address these questions in future work.

Remark 22. We note that the birth-death process studied above is a typical example of the type of setting where the N min -N max process may outperform the FV IPS. Here, we briefly outline another example where this is also the case.

We consider a branching random walk on {0, . . . , n} for some n ≥ 1 with a time inhomogeneous branching rate. More precisely, let (X t ) t≥0 denote a continuous time random walk on Z that moves one step to the right with probability p ∈ (0, 1) and one step to the left otherwise. Each step is taken after an (independent) exponentially distributed time with parameter 1. Now let b, k : {0, . . . , n} → [0, ∞) be bounded measurable functions and E b , E k be exponential random variables, that are independent of each other and X, and define

τ b 1 = inf{t > 0 : t 0 b(X s )ds > E b }, and
τ k 1 = inf{t > 0 : t 0 k(X s )ds > E k }. Further, set τ 1 = τ b 1 ∧ τ k 1 . If τ 1 = τ b 1
, a new particle is added to the system at position X τ b 1 , which continues the same behaviour as the original particle. If τ 1 = τ k 1 , the particle is removed from the system. We can then define τ 2 to be the first time after τ 1 that one of the particles alive at time τ 1 branches or is killed, and this process continues iteratively.

Letting N t denote the number of particles alive at time t ≥ 0, and {X i t : i = 1, . . . , N t } their positions, the branching random walk is defined as

Z t = Nt i=1 δ X i t , t ≥ 0.
We now consider a time-inhomogeneous version of the above process. We consider a sequence of switching times 0 = T off

0 < T on 1 < T off 1 < T on 2 < T off 2 < .
. . , defined via

T on i = T off i-1 + e i-1 (s on ), i ≥ 1,
and

T off i = T on i + e i (s off ), i ≥ 1,
where e i (s off ) and e i (s on ) are independent (of each other and everything else) exponential random variables with rate s off and s on , respectively. For T off i-1 ≤ t < T on i , Z t evolves as above but with b = 0, i.e. each particle moves as a random walk until it is (possibly) killed, until time T on i . At time T on i , we sample from an exponential random variable with rate B ∈ (0, ∞) and set b to be this value. For T on i ≤ t < T off i , Z t then evolves as the branching random walk with branching rate b for each particle. At time T off i , we again set the branching rate to be 0, and continue this process until some terminal time T .

We can see that in this case, while the branching rate is bounded almost surely on any compact time interval, we may simulate the process according to the Algorithm 1. However, we cannot use the fixed size Moran type algorithm since we do not know the upper bound for the branching rate.

However, there are other circumstances where the FV IPS may perform at least as well as the N min -N max process due to large deviations effects. We refer the reader to the next section for a discussion and possible solution to this.

Particle filtering methods for the penalised particle system

In many applications, it can be of interest to estimate numerically the growth rate of the semigroup Q for large times. A natural conjecture is that for large times T , m 0 Q T f ≈ m 0 , ψ ϕ, f e λT , where ϕ and ψ are the left eigenmeasure and right-eigenfunction respectively of the operator Q, and λ is the Lyapunov exponent for Q. One can try to estimate λ by computing the quantity m 0 Q T f for large times T .

By Theorem 6, one method for estimating λ would be through equation ( 6), so that we can simulate the (m T , Π A T , Π B T )-system and estimate λ by:

λ = 1 T log Π A T Π B T m T (1) = 1 T log Π A T Π B T , (33)
for T large. In practice, however, this is numerically a poor choice of estimator, due to largedeviation effects. The presence of the 'exponential' functionals, Π A T , Π B T mean that significant contributions to the expectation will be made by paths of the N min -N max process which have probability e -δT , and so infeasibly many simulations will be needed to compute an unbiased estimate for λ. A simple solution is to replace the estimator by

λ = 1 δt log 1 n i<n Π A t i+1 Π B t i+1 Π A t i Π B t i , (34)
where

t i = T i N and δt = T N .
In many examples, under conditions which ensure that the process is ergodic and thus converges to a stationary distribution, as T → ∞, for fixed δt, the estimator λ can be shown to converge to λ (e.g. [START_REF] Del Moral | Particle approximations of Lyapunov exponents connected to Schrödinger operators and Feynman-Kac semigroups[END_REF]).

In some cases the expression in [START_REF] Groisman | Simulation of quasi-stationary distributions on countable spaces[END_REF] may not be valid and we can only use estimates of the form in [START_REF] Grigorescu | Ergodic properties of multidimensional Brownian motion with rebirth[END_REF]. Examples may be when the process is not ergodic, or the system is not timehomogenous, such as the example given in Remark 22. In such circumstances, we propose a novel Monte Carlo method to compute estimates of λ. Specifically, we note that the N min -N max process X := {X 1 t , . . . , X Nt t } defines a continuous time Markov process, while the processes A t and B t corresponding to the number of selection and resampling events, are themselves additive functionals of the Markov process. Estimation of λ corresponds to computing a Feynman-Kac type exponential functional for a Markov process, and it is well known that estimating such functionals can be performed efficiently using interacting particle simulations, [START_REF] Del Moral | A Moran particle system approximation of Feynman-Kac formulae[END_REF][START_REF] Del Moral | Branching and interacting particle systems approximations of feynman-kac formulae with applications to non-linear filtering[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF][START_REF] Moral | Feynman-Kac formulae. Probability and its Applications[END_REF]. This suggests an algorithm as follows, assuming that b i , κ i , p i and q i have been chosen:

1. Fix T, N * , δt, N 0 , x 0 . Initiate with t = 0, and X 1 , . . . , X N * each containing N 0 particles of initial value x. Set W = 1.

2. Run each X i from time t to time t = min{t + δt, T } according to the N min -N max algorithm above, setting

w i := Π A i t Π B i t Π A i t Π B i t . 3. Set W := W × 1 N * N * i=1 w i .
4. If t < T , set t = t , resample each X i from the existing particles {X 1 , . . . , X N * } with probabilities proportional to w 1 , . . . , w N * , independently (multinomial resampling). Return to step 2.

If t = T , stop, and return

λ := 1 T log (W ) .
Here, the choice of multinomial resampling could be replaced by other resampling schemes. Similarly, the resampling could be conditional on properties of the particles (e.g. the effective sample size). We compare the effectiveness of the two methods in Figure 2, for a simple Markov process where the true value of the Lyapunov exponent can be calculated accurately numerically. We note that both methods have similar computational effort, and neither appears to converge to the 'true' eigenvalue, however the particle filter method achieves a much better numerical estimate of λ. The true value of λ is highlighted with the orange line, the green line represents an analytical estimate of the convergence value for the simulation (in this case, an estimate of the average growth of the exponential weights under the stationary measure for the system). In the second figure (right), λ is estimated using an interacting particle system, as described above. In this case, the particle filter has 100 particles, and is run for a total time of T = 40, giving comparable computational effort. It should be noted however that the second case also has a computational cost associated with the resampling scheme applied, and so the overall computational cost will be higher depending on the resampling scheme applied, and the frequency of resampling events.

Proof of Theorem 6

This section is dedicated to the proof of Theorem 6, which is rather long and so we break it up into five steps. In the first three steps, the main idea is to define two martingales M and M such that

ν f t -ν f 0 = Π A t Π B t (M t -M τ C t ) + At n=1 Π A ρn Π B ρn (M ρn -M ρn-) + Bt n=1 Π A σn Π B σn (M σn -M σn-) + Ct n=1 Π A τ n-1 Π B τ n-1 (M τn -M τ n-1 ), (35)
where σ n and ρ n denote the n th selection and resampling events, respectively, and

ν f t = Π A t Π B t i∈ St Q T -t f (X i t ).
The martingale M compensates the evolution of the branching process between resampling/selection times. The martingale M compensates the evolution of the particle system during resampling/selection events. The many-to-one formula (6) then follows. The final two steps make use of the following technical lemma to control the L 2 norms of the martingale increments in order to obtain [START_REF] Burdzy | A Fleming-Viot particle representation of the Dirichlet Laplacian[END_REF].

Lemma 23. For all ≥ 1, for all T ≥ 0,

E Π B T ≤ exp(c b ∞ T ), where c := sup x≥2 x (1 + 1/x) -1 .
In addition, denoting by β T the total number of branching events (with or without selection) up to time T , we have

E β T Π B T 2 ≤ N 0 exp((3 + c 4 ) b ∞ T /2).
The above estimates are the main reason why we require b to be bounded. These types of estimates are not difficult to derive in some particular situations where the branching rate is unbounded (typically for time inhomogeneous branching rates, or in situations where the branching rates is high in rarely visited regions), and in these situations, our proofs adapt easily. We defer the proof of the above lemma to the end of the section and now proceed to the proof of Theorem 6.

Proof of Theorem 6.

Step 1. The martingale M. For all i ∈ N, we define respectively the birth time B i and the death time D i of the particle with index i as

B i = inf{t ≥ 0, i ∈ St } and D i = inf{t ≥ B i , i / ∈ St }.
Fix T > 0 and t ∈ [0, T ]. For all particles i ∈ N, we set

M i,n t =      0, if t ≤ B i ∨ τ n , Q T -t f (X i t ) -Q T -B i ∨τn f (X i B i ∨τn ), if t ∈ [B i ∨ τ n , D i ∧ τ n+1 ), (1 -1 D n (i) + 1 B n (i) )Q T -D i ∧τ n+1 f (X i D i ∧τ n+1 ) -Q T -B i ∨τn f (X i B i ∨τn ), if t ≥ D i ∧ τ n+1 , (36)
where B n (i) denotes the event "τ n+1 is a branching time for the i-th particle" and D n (i) denotes the event "τ n+1 = D i is a soft killing time for the i-th particle".

By definition of the process between the time intervals [τ n , τ n+1 ), to show that M is indeed a martingale, it is sufficient to prove that

M 1,0 t := Q T -t f (X 1 t ) -Q T f (X 1 0 ), if t ∈ [0, τ 1 ), (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X τ 1 ) -Q T f (X 0 ), if t ≥ τ 1 ,
is a martingale. We observe that, by construction,

τ 1 = τ ∂ ∧ τ κ ∧ τ b where τ b = inf j∈ S0 inf{t ≥ 0, t 0 b j ((X i u ) i∈ Su ) du ≥ e b }, τ κ = inf j∈ S0 inf{t ≥ 0, t 0 κ j ((X i u ) i∈ Su ) ds ≥ e κ }, τ ∂ = inf j∈ S0 inf{t ≥ 0 X j t ∈ ∂}.
Then, for all s, t ≥ 0 such that s + t ≤ T , setting

E t = E[• | σ( Su , Xu , 0 ≤ u ≤ t)
] and using the fact that Q u f ≡ 0 on ∂ for all u ≥ 0, we have

E t (M 1,0 t+s ) = E t Q T -(t+s) f (X 1 t+s ) -Q T f (X 0 ) 1 t+s<τ 1 + E t x (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X τ 1 ) -Q T f (X 0 ) 1 t<τ 1 ≤t+s + (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X τ 1 ) -Q T f (X 0 ) 1 τ 1 ≤t (37)
Using the fact that the particles are independent copies of X before time τ 1 , and the Markov property at time t, we obtain

E t x Q T -(t+s) f (X 1 t+s )1 t+s<τ 1 = E t Q T -(t+s) f (Y 1 t+s )e -t+s t hu du 1 t+s<τ ∂ 1 t<τ 1 ,
where E t is the expectation on a probability space Ω such that (Y i u ) u≥t , i ∈ S0 , are independent copies of X, starting from X i t at time t, and where

h u := j∈ S0 b j ((Y i u ) i∈ S0 ) + κ j ((Y i u ) i∈ S0 ).
Then, using the Markov property at time t + s, we obtain

E t x Q T -(t+s) f (X 1 t+s )1 t+s<τ 1 = E t e T t+s b(Y 1 u )-κ(Y 1 u ) du f (Y 1 T )e -t+s t hu du 1 t+s<τ ∂ 1 t<τ 1 = E t e T t b(Y 1 u )-κ(Y 1 u ) du f (Y 1 T )e -t+s t h 1 u du 1 t+s<τ ∂ 1 t<τ 1 ,
where

h 1 u := 2b 1 ((Y i u ) i∈ S0 )+ j∈ S0 ,j =1 b j ((Y i u ) i∈ S0 )+κ j ((Y i u ) i∈ S0
). Similarly, since κ 1 is the intensity of the soft killing time for the first particle, and b 1 is the intensity of its branching time, we obtain

E t x (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X 1 τ 1 )1 t<τ 1 ≤t+s = E t (t+s)∧τ ∂ t (h u -κ 1 ((Y i u ) i∈ S0 ) + b 1 ((Y i u ) i∈ S0 ))Q T -u f (Y 1 u )e -u t hv dv du 1 t<τ 1 + E t (1 -1 Y 1 τ ∂ ∈∂ )e -τ ∂ t hu du 1 τ ∂ ≤t+s Q T -τ ∂ f (Y 1 τ ∂ ) 1 t<τ 1 .
Hence, using the fact that Y satisfies the strong Markov property at time τ ∂ , and the fact that

Q T -τ ∂ f (Y 1 τ ∂ ) = 0 on the event {Y 1 τ ∂ ∈ ∂}, we deduce that E t x (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X 1 τ 1 )1 t<τ 1 ≤t+s = E t (t+s)∧τ ∂ t h 1 u e T u (b(Y 1 v )-κ(Y 1 v )) dv f (Y 1 T )e -u t hv dv du 1 t<τ 1 + E t e -τ ∂ t hu du 1 τ ∂ ≤t+s e T τ ∂ (b(Y 1 v )-κ(Y 1 v )) dv f (Y 1 T ) 1 t<τ 1 = E t e T t (b(Y 1 v )-κ(Y 1 v )) dv f (Y 1 T ) (t+s)∧τ ∂ t h 1 u e -u t h 1 v dv du + e -τ ∂ t h 1 v dv 1 τ ∂ ≤t+s 1 t<τ 1 ,
where we used the assumption b(Y

1 v ) -κ(Y 1 v ) = b 1 ((Y i u ) i∈ S0 ) -κ 1 ((Y i u ) i∈ S0 ). But e -t+s t h 1 u du 1 t+s<τ ∂ + (t+s)∧τ ∂ t h 1 u e -u t h 1 v dv du + e -τ ∂ t h 1 v dv 1 τ ∂ ≤t+s = 1
Hence we deduce from (37) that

E t x (M 1,0 t+s ) = (1 -1 D(1) + 1 B(1) ) (Q T -τ 1 f (X τ 1 ) -Q T f (X 0 )) 1 τ 1 ≤t + E t e T t (b(Y 1 v )-κ(Y 1 v )) dv f (Y 1 T ) -Q T f (x) 1 t<τ 1 = (1 -1 D(1) + 1 B(1) )Q T -τ 1 f (X τ 1 )1 τ 1 ≤t + Q T -t f (X 1 t ) 1 t<τ 1 -Q T f (x) = M 1,0 t .
This concludes Step 1.

Step 2. A first martingale decomposition. Fix t ≥ 0 and recall that C t = inf{n ≥ 0 : t < τ n+1 }. In what follows, we denote by ρ 1 , ρ 2 , . . . the sequence of times at which a resampling event occurs, and by σ 1 , σ 2 , . . . the sequence of times at which a selection event occurs. Note that, if τ n ∈ {ρ 1 , ρ 2 , . . .} or τ n ∈ {σ 1 , σ 2 , . . .}, then N n = N n-1 ∈ {N min , . . . , N max }.

We will show that, for all t ≥ 0, we have

M t := Ct n=0 i∈ St M i,n t = i∈ St Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ) - At n=1 Q T -ρn f (X in ρn ) + Bt n=1 Q T -σn f (X j n-1 σn ), (38)
where we recall that A t is the number of resampling events up to time t, B t is the number of selection events up to time t, and where we have set i n to denote the index of the particle added at time ρ n and j n-1 is the particle killed (due to a selection event) at time σ n (as above, X j n-1 ρn denotes the position of the particle before the selection event). We will show that [START_REF] Ikeda | Branching Markov processes II[END_REF] holds for n = 0, 1, since the general formula then follows using similar arguments. Noting that for all t ∈ [τ 0 , τ 1 ), C t = A t = B t = 0, in this case we have

Ct n=0 i∈ St M i,n t = i∈ S0 M i,0 t = i∈ S0 Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ),
as required.

Let us now describe the evolution of the martingale at time τ 1 , depending on the type of event that occurs.

Case τ 1 = ρ 1 . In this case, for all t ∈ [τ 1 , τ 2 ), we have A t = 1 and B t = 0, and St = {1 + max S0 } ∪ S0 \ {i 0 }, where i 0 is the particle removed at time τ 1 . Thus

Ct n=0 i∈ St M i,n t = i∈ S0 M i,0 t + i∈ Sτ 1 M i,1 t = i∈ S0 (1 -1 i=i 0 )Q T -τ 1 f (X i τ 1 ) -Q T f (X i 0 ) + i∈ Sτ 1 Q T -t f (X i t ) -Q T -τ 1 f (X i τ 1 ) = i∈ Sτ 1 Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ) -Q T -τ 1 f (X 1+max S0 τ 1 ) = i∈ St Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ) -Q T -τ 1 f (X i 1 τ 1 ),
as claimed.

Case τ 1 = σ 1 . Similarly, if τ 1 = σ 1 , then, for all t ∈ [τ 1 , τ 2 ), we have A t = 0 and B t = 1, and St = {1 + max S0 } ∪ S0 \ {j 0 }. Denoting by j 0 the particle duplicated at time τ 1 , we have

Ct n=0 i∈ St M i,n t = i∈ S0 M i,0 t + i∈ Sτ 1 M i,1 t = i∈ S0 (1 + 1 i=j 0 )Q T -τ 1 f (X i τ 1 ) -Q T f (X i 0 ) + i∈ Sτ 1 Q T -t f (X i t ) -Q T -τ 1 f (X i τ 1 ) = i∈S 1 Q T -t f (X i t ) - i∈S 0 Q T f (X i 0 ) + Q T -τ 1 f (X j 0 τ 1 ) = i∈ St Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ) + Q T -τ 1 f (X j 0 τ 1 ).
Case τ 1 / ∈ {ρ 1 , σ 1 }. If τ 1 / ∈ {ρ 1 , σ 1 }, we have, for all t ∈ [τ 1 , τ 2 ), A t = B t = 0. If τ 1 is a branching event of the particle with index j 0 ∈ S 0 , then St = {1 + max S0 } ∪ S0

Ct n=0 i∈ St M i,n t = i∈ S0 (1 + 1 i=j 0 )Q T -τ 1 f (X i τ 1 ) -Q T f (X i 0 ) + i∈ Sτ 1 Q T -t f (X i t ) -Q T -τ 1 f (X i τ 1 ) = i∈ Sτ 1 Q T -t f (X i t ) - i∈ S0 Q T f (X i 0 ),
where we used the fact that

X j 0 τ 1 = X 1+max S0 τ 1
. Continuing in this manner, we obtain [START_REF] Ikeda | Branching Markov processes II[END_REF].

Step 3. A second martingale decomposition.

For all time t ≥ 0 and all bounded measurable function f : E → R, we denote the integral of f with respect to the occupation measure of the particle system at time t by

µ f t := i∈ St Q T -t f (X i t ). ( 39
)
Further, recall the definitions of Π A t , Π B t and

ν f t = Π A t Π B t µ f t .
In this step, we prove that equation [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF] holds with M defined as in the previous steps, and M to be defined shortly. To this end, we first note that from Step 2 we obtain

µ f t -µ f 0 = Ct n=0 i∈ St M i,n t + At n=1 Q T -ρn f (X i n-1 ρn ) - Bt n=1 Q T -σn f (X j n-1 σn ). (40)
For n ∈ N, consider the following martingale increments. We define

M ρn -M ρn-:= Q T -ρn f (X i n-1 ρn ) - 1 Nρn -1 i∈ Sρn-\{i n-1 } Q T -ρn f (X i ρn ), (41)
where i n-1 is the index of the particle duplicated at time ρ n and i n-1 is the index of the particle removed at time ρ n (note also that Nρn = Nρn-= | Sρn-|). We also define

M σn -M σn-:= 1 Nσn + 1   i∈ Sσn- Q T -σn f (X i σn ) + Q T -σn f (X j n-1 σn )   -Q T -σn f (X j n-1 σn ), ( 42
)
where we recall that j n-1 is the particle removed at time σ n and j n-1 is the index of the particle duplicated during the branching event at this time (note also that Nσn = Nσn-= | Sσn-|). We assume further that M is constant except at these jump times. Then,

µ f t -µ f 0 = M t + M t + At n=1 1 Nρn -1 i∈ Sρn-\{i n-1 } Q T -ρn f (X i ρn ) - Bt n=1 1 Nσn + 1   i∈ Sσn- Q T -σn f (X i σn ) + Q T -σn f (X j n-1 σn )   . ( 43
)
For each n ≥ 1, we define µ f τn-and ν f τn-by

µ f τn-= i∈S n-1 f (X i τn ) -1 τn∈K f (X i n-1 τn ) + 1 τn∈B f (X j n-1 τn ) and ν τn-= Π A τ n-1 Π B τ n-1 µ τn-
where K is the set of killing times and B is the set of branching times. Informally, µ τn-and ν τn-represent the state of the particle system at time τ n before the resampling or selection eventually occurring at time τ n . Note that, at any time τ n that is not a resampling or a selection time, we have µ τn-= µ τn and ν τn-= ν τn . Hence, denoting by (θ n ) n the sequence of events from the families ρ or σ, we obtain for all t ≥ 0,

ν f t -ν f 0 = ν t -ν τ C t + Ct n=1 (ν τn -ν τn-) + Ct n=1 ν τn--ν τ n-1 = Π A t Π B t (µ t -µ τ C t ) + At n=1 ν ρn -ν ρn-+ Bt n=1 ν σn -ν σn-+ Ct n=1 ν τn--ν τ n-1 .
At resampling times ρ n , we have Π B ρn -Π B ρn-= 0 and hence,

ν f ρn -ν f ρn-= Π A ρn Π B ρn µ f ρn -Π A ρ n-1 Π B ρ n-1 µ f ρn- = Π A ρn Π B ρn (µ f ρn -µ f ρn-) + Π B ρn µ f ρn-Π A ρn -Π A ρ n-1 = Π A ρn Π B ρn Q T -ρn f (X i n-1 ρn ) + Π B ρn µ f ρn-Π A ρn -Π A ρ n-1 (44)
For the increment of Π A , we note that

Π A ρn -Π A ρ n-1 = Π A ρn - Nρn Nρn -1 Π A ρn = - 1 Nρn -1 Π A ρn .
Substituting this into (44) and using [START_REF] Jagers | Branching processes as population dynamics[END_REF], we obtain

ν f ρn -ν f ρn-= Π A ρn Π B ρn Q T -ρn f (X i n-1 ρn ) - 1 Nρn -1 µ f ρn- = Π A ρn Π B ρn (M ρn -M ρn-) . (45)
Similarly, for n ∈ {1, . . . , B t },

ν f σn -ν f σn-= Π A σn Π B σn (µ f σn -µ f σn-) + Π A σn µ f σn-Π B σn -Π B σn-. (46)
In this case, we have

Π B σn -Π B σn-= Π B σn - Nσn Nσn + 1 Π B σn = 1 Nσn + 1 Π B σn , (47)
and

µ f σn -µ f σn-= -Q T -σn f (X j n-1 σn ), (48)
where j n-1 is the index of the particle removed from the system at time σ n . Again, substituting these equalities back into [START_REF] Li | Branching processes with immigration and related topics[END_REF] and using [START_REF] Kawazu | Branching processes with immigration and related limit theorems[END_REF] gives,

ν f σn -ν f σn-= Π A σn Π B σn -Q T -σn f (X j n-1 σn ) + 1 Nσn + 1 µ f σn- = Π A σn Π B σn (M σn -M σn-). (49)
Finally, we have ν τn--

ν τ n-1 = Π A τ n-1 Π B τ n-1 (M τn -M τ n-1
) for all n ≥ 1, hence

ν f t -ν f 0 = Π A t Π B t (M t -M τ C t ) + At n=1 Π A ρn Π B ρn (M ρn -M ρn-) + Bt n=1 Π A σn Π B σn (M σn -M σn-) + Ct n=1 Π A τ n-1 Π B τ n-1 (M τn -M τ n-1 ),
which entails [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF] and thus concludes Step 3.

Step 4. Control of the L 2 norm of the martingale increments. In this section, we assume without loss of generality (taking

e -b ∞T f / f ∞ instead of f ) that sup t∈[0,T ] Q t f ∞ ≤ 1.
We now show that the L 2 norm of each of these quantities on the righthand side of ( 35) is finite. To this end, let us first consider M ρn -M ρn-. We have, according to [START_REF] Jagers | Branching processes as population dynamics[END_REF],

|M ρn -M ρn-| = Q T -ρn f (X i n-1 ρn ) - 1 Nρn -1 i∈ Sρn-\{i n-1 } Q T -ρn f (X i ρn ) ≤ 2. (50)
Then, setting

R T = max t∈[0,T ] N t , we obtain E   A T n=1 Π A ρn Π B ρn (M ρn -M ρn-) 2   = E   Nρn -1 Nρn A T n=1 Π A ρn-Π B ρn-(M ρn -M ρn-) 2   ≤ 4 E A T n=1 R T -1 R T 2Aρ n (Π B ρn-) 2 ≤ 4 E (Π B T ) 2 A T n=1 R T -1 R T 2n . (51)
To control the sum term, since 1/R T ≤ 1 2 , we have

A T n=1 1 - 1 R T 2n ≤ +∞ n=1 1 - 1 R T 2n = (1 -1 R T ) 2 1 -(1 -1 R T ) 2 = R T 1 -2/R T + R -2 T 2 -1 R T ≤ R T /2 ≤ (N 0 + β T )/2,
where β T denotes the total number of branching events before time T . We now make use of Lemma 23, stated and proved at the end of this section, which entails that

E (Π B T ) 2 (N 0 + β T )/2 ≤ N 0 exp(c b ∞ T ), (52)
for some (explicit) constant c > 0. Substituting these estimates into (51), we obtain

E   A T n=1 Π A ρn Π B ρn (M ρn -M ρn-) 2   ≤ 4 N 0 exp (c b ∞ T ) . (53)
Let us now consider ∆ σn M. Using (42), we have

|M σn -M σ n-| ≤ 2. Then E   B T n=1 Π A σn Π B σn (M σn -M σ n-) 2   (54) ≤ E   B T n=1 Π A σ n-Π B σ n- N σ n-+ 1 N σ n- (M σn -M σ n-) 2   ≤ 4 E B T n=1 Π B σn 2 ≤ 4 E β T Π B T 2 . ( 55
)
Using Lemma 23 stated and proved at the end of this section, we deduce that

E   B T n=1 Π A σn Π B σn (M σn -M σ n-) 2   ≤ 4 N 0 exp(c b ∞ T ), (56)
where c can be chosen to be the same as in [START_REF] Tough | The Fleming-Viot process with McKean-Vlasov dynamics[END_REF].

It remains to control the increments of M. For all n ≥ 1, let E n ⊂ N be the set of indices of particles defined by E n = {i ∈ N s.t. the i th particle is branching, killed, born, removed or duplicated at time τ n }.

Observe that |E n | ≤ 4 for all n ≥ 1. We also define, for all i ∈ N, the sequence of times (τ i n ) n≥0 by τ i 0 = 0 and

τ i n+1 = min{τ k > τ i n , with k ≥ 0 and i ∈ E k }.
We write

M T = i∈N n≥1 M i τ i n ∧T -M i τ i n-1 ∧T .
Since all the martingales increments are orthogonal, we have, for all i ∈ N,

∆ := E   Π A T Π B T (M T -M τ C T ) + C T n=1 Π A τ n-1 Π B τ n-1 (M τn -M τ n-1 ) 2   = E   i∈N n≥1 Π A τ i n-1 Π B τ i n-1 (M i τ i n ∧T -M i τ i n-1 ∧T ) 2   = E i∈N n≥1 Π A τ i n-1 Π B τ i n-1 2 M i τ i n ∧T -M i τ i n-1 ∧T 2 .
But the increments M i τ i n+1 ∧T -M i τ i n ∧T are almost surely bounded by 3, so that

∆ ≤ 9 E i∈N n≥1 Π A τ i n-1 Π B τ i n-1 2 1 τ i n-1 <T ≤ 9 E n≥1 i∈En Π A τ n-1 Π B τ n-1 2 1 τ n-1 <T ≤ 36 E n≥1 Π A τ n-1 Π B τ n-1 2 1 τ n-1 <T
Hence, denoting by β T (resp. κ T ) the number of branching (resp. killing without resampling) events occuring before time T ≥ 0 (so that

C T = β T + κ T + A T ), we have ∆ ≤ 36 E Π B T 2 (1 + β T + κ T + A T n=1 Π A ρ n-1 2 ) ( 57
)
We already computed that, almost surely,

A T n=1 Π A ρ n-1 2 ≤ +∞ k=0 R T -1 R T 2k ≤ R T /2 ≤ (N 0 + β T )/2.
Now, observe that the number of particles in the system as time T is upper bounded by N 0 + β T -κ T and lower bounded by 0, so that κ T ≤ N 0 + β T , and hence

E Π B T 2 (1 + β T + κ T + A T n=1 Π A ρ n-1 2 ) ≤ E Π B T 2 (β T + 3N 0 ) ≤ 4N 0 exp(c b ∞ T ),
where we used, as above, Lemma 23. We thus obtained, using the last inequality and [START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF], that

∆ ≤ 144 N 0 exp (c b ∞ T ) .
Using [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF], ( 53), [START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF] and the last inequality, we deduce that

E ν f T -ν f 0 2 ≤ 152N 0 exp (c b ∞ T ) .
Without assuming that sup

t∈[0,T ] Q t f ∞ ≤ 1, we thus obtain E ν f T -ν f 0 2 ≤ 152N 0 exp (c b ∞ T ) sup t∈[0,T ] Q t f ∞ 2 . ( 58
)
Step 5. Conclusion.

Recall the occupation measure m t defined in (5) for all t ≥ 0. With this notation, the conclusion of Step 4 reads

E m 0 Q T f -Π A T Π B T m T f 2 ≤ 152N 0 exp (c b ∞ T ) sup t∈[0,T ] Q t f ∞ 2 .
We deduce that

m 0 Q T 1 E m T f m T 1 E 1 m T =0 -m 0 Q T f 2 ≤ (m 0 Q T 1 E -Π A T Π B T m T 1 E ) m T f m T 1 E 1 m T =0 2 + Π A T Π B T m T f -m 0 Q T f 2 ≤ 2 152N 0 exp (c b ∞ T /2) f ∞ sup t∈[0,T ] Q t 1 E ∞ .
We conclude that

m 0 Q T f m 0 Q T 1 E - m T f m T (1 E ) 1 m T =0 2 ≤ 2 √ 152 exp (c b ∞ T /2) f ∞ sup t∈[0,T ] Q t 1 E ∞ m 0 Q T 1 E /N 0 1 √ N 0 .
Observing that sup

t∈[0,T ] Q t 1 E ∞ ≤ e b ∞T ,
this concludes the proof of Theorem 6.

Proof of Lemma 23. We use the same notations as in the proof of Theorem 6.

Fix C > c b ∞ and δ > 0 such that C -δC 2 2 ≥ c b ∞ . We have E Π B τ 1 ∧δ e -C(τ 1 ∧δ) = E N 0 + 1 N 0 1 σ 1 =τ 1 ∧δ e -C(τ 1 ∧δ) + E 1 σ 1 >τ 1 ∧δ e -C(τ 1 ∧δ) = N 0 + 1 N 0 E e -C(τ 1 ∧δ) + 1 - N 0 + 1 N 0 E 1 σ 1 >τ 1 ∧δ e -C(τ 1 ∧δ) .
On the one hand, we have

E e -C(τ 1 ∧δ) ≤ 1 -e -Cδ CE(τ 1 ∧ δ) ≤ 1 -(C -δC 2 )E(τ 1 ∧ δ),
and, on the other hand,

E 1 σ 1 >τ 1 ∧δ e -C(τ 1 ∧δ) ≥ E e -(τ 1 ∧δ)(N 0 b ∞+C) ≥ 1 -(N 0 b ∞ + C)E (τ 1 ∧ δ) .
Setting, for the sake of readability, α = E(τ 1 ∧ δ) and β = N 0 + 1 N 0 , we deduce that

E Π B τ 1 ∧δ e -C(τ 1 ∧δ) ≤ 1 -α C -δC 2 β + N 0 b ∞ (1 -β) ≤ 1 -α C -δC 2 2 -b ∞ c ≤ 1, (59)
where we used the fact that β ≤ 2 and β -1 ≤ c /N 0 . We define the sequence of stopping times (θ n ) n≥0 by θ 0 = 0 and

θ n+1 := inf{((n + 1)δ) ∧ τ k , k ≥ 1, τ k > θ n }. Then E Π B T ∧θ n+1 e -C(T ∧θ n+1 ) | F T ∧θn = Π B T ∧θn e -C(T ∧θn) × E X T ∧θn Π B (T -u)∧θ 1 e -C((T -u)∧θ 1 )
|u=T ∧θn .

Since T -u ≤ δ almost surely in the right hand term, we deduce from ( 59) that

E Π B T ∧θ n+1 e -C(T ∧θ n+1 ) | F T ∧θn ≤ Π B
T ∧θn e -C(T ∧θn) .

In particular, using Fatou's Lemma, we deduce, letting n → +∞, that

E Π B T e -CT ≤ 1.
Since this is true for any C > c b ∞ , this concludes the proof of the first part of the lemma.

For the second part, we observe that β T is stochastically dominated by a continuous time counting process (Z t ) t≥0 on N starting from N 0 and with increment rate b ∞ n from state n to state n + 1. In particular, denoting by L the infinitesimal generator of this process and V (n) = n 2 for all n ≥ 1, we obtain

LV (n) = b ∞ n((n + 1) 2 -n 2 ) = b ∞ n(2n + 1) ≤ 3 b ∞ V (n), and deduce that E(β 2 T ) ≤ E(Z 2 T ) ≤ e 3 b ∞T Z 2 0 = e 3 b ∞T N 2 0 . Hence, obtain E β T Π B T 2 ≤ E(β 2 T ) E (Π B T ) 4 ≤ N 0 exp((3 + c 4 ) b ∞ T /2).
This concludes the proof of the second part of the Lemma.

Remark 24. Our proof adapts to the case where τ ∞ < +∞ with positive probability (that is when Assumption 2 does not holds true), using the following strategy. In this situation, one shows using the same approach that, for all n ≥ 0,

m 0 Q T f = E Π A τn∧T Π B τn∧T m τn∧T Q T -τn∧T f .
Since the branching rate is uniformly bounded and since Π A τn∧T is bounded by 1, one checks that Π A τn∧T Π B τn∧T m τn∧T Q T -τn∧T f can be uniformly bounded (in n) by an integrable random variable. Since in addition Π A τn∧T Π B τn∧T m τn∧T Q T -τn∧T f → 0 on the event τ ∞ ≤ T (since this events necessarily corresponds to the accumulation of resampling events) and

Π A τn∧T Π B τn∧T m τn∧T Q T -τn∧T f → Π A T Π B
T m T f on the complementary event, one deduces from the dominated convergence theorem that

m 0 Q T f = E Π A T Π B T m T f 1 T <τ∞ .
The rest of the proof (see Steps 4 and 5) then proceeds almost identically.

Remark 25. We comment on some direct and natural continuations of the proof, which are easily derived from the fact that the proof is based on the control of martingale increments. First, if the martingale ν f t -ν f 0 is càdlàg, then Doob's maximal inequality and [START_REF] Yanev | Critical controlled branching processes and their relatives[END_REF] entail that

E sup t∈[0,T ] ν f T -ν f 0 2 ≤ c T f 2 ∞ sup t∈[0,T ] Q t 1 E 2 ∞ N 0
for some (explicit) constant c T . Then the conclusion reads

sup t∈[0,T ] m 0 Q t f m 0 Q t 1 E - m t f m t (1 E ) 1 mt =0 2 ≤ c T f ∞ m 0 Q T 1 E /N 0 1 √ N 0 ,
Second, if the killing rate is bounded, then the statement of Lemma 23 holds true with β T replaced by the total number of killing/branching events C T , and with different constants. Then, for any p ≥ 2, using the same approach as in Step 4 of the proof, but using the Burkholder-Davis-Gundy inequality instead, one deduces similarly that

E ν f T -ν f 0 p ≤ c T N p/2 0 ,
and hence

m 0 Q t f m 0 Q t 1 E - m t f m t (1 E ) p ≤ c T f ∞ m 0 Q T 1 E /N 0 1 √ N 0 ,
for some (explicit) constants c T and c T . Note that the assumption that the killing rate is uniformly bounded is technical: after using Burkholder-Davis-Gundy inequality, one needs to control the term

E C p/2 T N + 1 N pC T
, however, depending on the particular process at hand, there may be other ways to obtain upper bounds on this quantity which do not require the killing rate to be bounded.

A Formal construction of the particle system

In this section, we provide a formal construction of the particle system described above, and show that it defines a Markov process with respect to its natural filtration. We consider the measurable space

F = s∈P f (N) {s} ⊗ Ω s ⊗ R s + ⊗ R s + ,
where, as usual, denumerable spaces are endowed with the discrete σ-fields, R + is endowed with its Borel σ-field, and product spaces are endowed with the product σ-field. We denote generic elements of F by θ = (s, (ω i ) i∈s , (e b i ) i∈s , (e κ i ) i∈s ), where s ∈ P f (N) is used to enumerate the particle system, where, for all i ∈ s, X(ω i ) : t ∈ [0, +∞) → X t (ω i ) ∈ E ∪ ∂ is the trajectory of the i th particle, and where e b i ∈ R + and e κ i ∈ R + are auxiliary components used to define the branching and soft-killing times.

Fix s ∈ P f (N) and let us consider the probability kernel N 1 from {s} ⊗ E s to F defined by N 1 (s, (x i ) i∈s ) := δ s ⊗ i∈s

P x i ⊗ E ⊗s ⊗ E ⊗s ,
where E denotes the exponential distribution with parameter 1. This kernel is used to model the evolution of the particles before the occurence of a branching/killing event.

Given an element θ ∈ F , our aim is now to define the position of the system at the first time of branching/killing. In order to do so, we define, for all θ ∈ F and all i ∈ s, the functions τ i,b (θ) = inf{t ≥ 0, t 0 b i 0 ((X u (ω j )) j∈s ) du = e b i }, τ i,κ (θ) = inf{t ≥ 0, t 0 κ i 0 ((X u (ω j )) j∈s ) du = e κ i }, τ i,∂ (θ) = inf{t ≥ 0, X t (ω i ) ∈ ∂}.

Note that, for all t ≥ 0, (u, ω) ∈ [0, +∞)×F → X u (ω) is measurable with respect B([0, +∞))⊗ F, which entails that ω → t 0 b i ((X u (ω j )) j∈s ) du is F ⊗s -measurable. In addition, since t → t 0 b i ((X u (ω j )) j∈s ) du, t → t 0 κ i ((X u (ω j )) j∈s ) du and t → 1 Xt(ω i )∈∂ are right-continuous, the infima in the above definitions can be taken over Q and are thus measurable functions of (ω j ) j∈s , (e b j ) j∈s and (e κ j ) j∈s and thus of θ. This is also the case for the functions τ and i 0 defined as τ (θ) := min i∈s τ i,b (θ) ∧ τ i,κ (θ) ∧ τ i,∂ (θ), i 0 (θ) := min{i ∈ S, τ i,b (θ) ∧ τ i,κ (θ) ∧ τ i,∂ (θ) = τ (θ)}.

For all θ ∈ F , we define the final position of θ, denoted by Y (θ) ∈ (E ∪ ∂) s , by Y (θ) := (X τ (ω) (ω i )) i∈s , which is also a measurable function of θ.

Then we construct a kernel which, to any given θ ∈ F and its final position Y (θ), associates its position after the branching/killing and selection/resampling event. In order to do so, we define the kernel N 2 from F to ∪ s∈P f (N) {s} × (E ∪ ∂) s by N 2 (θ) = 1 τ (θ)=τ i 0 ,b p b,i 0 (θ) (Y (θ)) 1 |s| + 1 i∈s + δ s -i + ⊗ δ Y (θ) +i 0 ,-i + 1 τ (θ)=τ i 0 ,b 1 -p b,i 0 (θ) (Y (θ)) δ s + ⊗ δ Y (θ) +i 0 + 1 τ (θ)=τ i 0 ,κ ∧τ i 0 ,∂ p κ,i 0 (θ) (Y (θ)) 1 |s| -1 i∈s -i 0 δ s -i 0 + ⊗ δ Y (θ) +i,-i 0 + 1 τ (θ)=τ i 0 ,κ ∧τ i 0 ,∂ 1 -p κ,i 0 (θ) (Y (θ)) δ s -i 0 ⊗ δ Y (θ) -i 0 where s + = {1 + max s}, s -i + = {1 + max s} ∪ s \ {i}, Y (θ) i 0 ,-i ∈ (E ∪ ∂) s -i + is such that Y (θ) i 0 ,-i j = Y (θ) j for all j ∈ s \ {i}, Y (θ) i 0 for j = 1 + max s, Y (θ) i 0 ∈ (E ∪ ∂) s + is such that Y (θ) i 0 j = Y (θ) j for all j ∈ s, Y (θ) i 0 for j = 1 + max s, Case m = n. In the former case, we observe that, by definition of N 1 , f (a, u) = P a X s +1 -u ∈ A +1 , . . . , X s k -u ∈ A k , s i -u ∈ [0, τ 1 ), X t-u ∈ A, t -u ∈ [0, τ 1 ) = P a Y s +1 -u ∈ A +1 , . . . , Y s k -u ∈ A k , s i -u ∈ [0, t), Y t-u ∈ A, t -u ∈ [0, t) , where P a is the product measure Ω = P ⊗s ⊗ E(1) s on the probability space Ω s (R + ) s ⊗ (R + ) s endowed with the filtration defined by H t := F ⊗s t ⊗ {∅, R + } s ⊗ {∅, R + } s , with the variables Y and t being defined, for all θ = ((ω i ) i∈s , (e b i ) i∈s , (e κ i ) i∈s ) ∈ Ω , by Y t = (s, (X t (ω i ) i∈s ), ∀t ≥ 0, and t = min i∈s t b,i ∧ t κ,i ∧ t ∂,i , with t b,i = inf{t ≥ 0, t 0 b i 0 ((X u (ω j )) j∈s ) du = e b i }, t κ,i = inf{t ≥ 0, t 0 κ i 0 ((X u (ω j )) j∈s ) du = e κ i }, t ∂,i = inf{t ≥ 0, X t (ω i ) ∈ ∂}.

We emphasise that the second component of Y is distributed as a vector indexed by s of independent copies of X starting from x i , for all i ∈ s, and that it is a progressively measurable Markov process with respect to the filtration H. In addition, the random variables e b i and e κ i , i ∈ s, are independent from each others and from H t , for all t ≥ 0.

Using the Markov property at time s -u for the process Y, we obtain where h(s, (x i ) i∈s ) := i∈s κ i (x i ) + b i (x i ), and where we used the fact that Y is progressively measurable to deduce that s-u 0 h(Y v ) dv is H s-u -measurable. We deduce that f (a, u) = E a 1 X s +1 -u ∈A +1 ,...,X s k -u ∈A k , s i -u∈[0,τ 1 ) E X s-u (X t-s ∈ A, t -s ∈ [0, τ 1 ))

P a Y s +1 -u ∈ A +1 , . . . , Y s k -u ∈ A k , s i -u ∈ [0, t), Y t-u ∈ A, t -u ∈ [0, t) = E a 1 Y s +1 -
Taking the expectation in (61), we finally deduce that Case m < n. We can now proceed with the case m < n. In this case, we write f (a, u) = E a 1 X s +1 -u ∈A +1 ,...,X s k -u ∈A k , s i -u∈[0,τ 1 ) g u (A 1 , τ 1 ) , where g u (b, v)

P µ 0 (X s 1 ∈ A 1 , . . . , X s k ∈ A k , s i ∈ [τ n i -1 , τ n i ), X s+t ∈ A, s + t ∈ [τ n , τ n+1 )) = E µ 0 1
:= E A 1 (X t-u-v ∈ A, t -u -v ∈ [τ n-m-1 , τ n-m ))
Using the same strategy as above (with the same definition of Y), we obtain, using also the Markov property at time s k -u = s -u, f (a, u) = E a 1 Y s +1 -u ∈A +1 ,...,Y s k -u ∈A k , s i -u∈[0,t) g u (Y t , t)

= E a 1 Y s +1 -u ∈A +1 ,...,Y s k -u ∈A k , s i -u∈[0,t) E Y s-u (g u (Y t , t))

= E a 1 X s +1 -u ∈A +1 ,...,X s k -u ∈A k , s i -u∈[0,τ 1 ) E X s-u (g u (X τ 1 , τ 1 ))

Using again the definition of the law of (A 0 , Θ Conclusion. Finally, using (62) and ( 63) and summing over the indices n i and n as in [START_REF] Yanev | Branching processes in cell proliferation kinetics[END_REF], we deduce that

P µ 0 (X s 1 ∈ A 1 , . . . , X s k ∈ A k , X t ∈ A) = E µ 0 1 Xs 1 ∈A 1 ,...,Xs k ∈A k P Xs (X t-s ∈ A) .
Since the property holds true for all finite family of time indices s 1 < . . . < s k = s and all A 1 , . . . , A k , we deduce that P µ 0 (X t ∈ A | σ(X u , u ≤ s)) = P Xs (X t-s ∈ A) , which concludes the proof of the proposition.

Remark 28. Although we won't use this property, the process ( St , (X i t ) t∈ St ) t≥0 can be shown to be a Markov process with respect to larger filtrations. On a related note, it is natural to ask whether the strong Markov property for X entails the strong Markov property for the particle system. This is not clear at all in full generality and we leave it as an open problem. However, it can be shown, using a similar approach as in [START_REF] Ikeda | Branching Markov processes III[END_REF], that, under additional regularity properties on the process X and on the functionals b, κ, and additional requirements on the state space E, the particle system is indeed a strong Markov process in a quite general setting.

Figure 1 :

 1 Figure 1: Simulations of the BBMMI model where X is a Brownian motion, killed at 0 and 1, with constant branching rate.

Figure 2 :

 2 Figure2: We use two methods described above to estimate the eigenvalue λ. The first simulation (left) estimates λ using[START_REF] Grigorescu | Ergodic properties of multidimensional Brownian motion with rebirth[END_REF], based on a single trajectory run up to time T = 4000. The true value of λ is highlighted with the orange line, the green line represents an analytical estimate of the convergence value for the simulation (in this case, an estimate of the average growth of the exponential weights under the stationary measure for the system). In the second figure (right), λ is estimated using an interacting particle system, as described above. In this case, the particle filter has 100 particles, and is run for a total time of T = 40, giving comparable computational effort. It should be noted however that the second case also has a computational cost associated with the resampling scheme applied, and so the overall computational cost will be higher depending on the resampling scheme applied, and the frequency of resampling events.

Table 2 :

 2

	10	N min -N max	0.01	0.12	144
		FV	0.02	0.18	857
	M = 100	N min -N max	0.01	0.12	144
		FV	0.10	0.39	9866
	M = 1000	N min -N max	0.01	0.12	144
		FV	0.20	0.50	99873
	M = +∞	N min -N max	0.01	0.12	144
		FV	*	*	*

  u ∈A +1 ,...,Y s k -u ∈A k 1 Y t-u ∈A\∂ e -t-u 0 h(Yv) dv = E a 1 Y s +1 -u ∈A +1 ,...,Y s k -u ∈A k E a 1 Y t-u ∈A\∂ e -t-u Y s +1 -u ∈A +1 ,...,Y s k -u ∈A k 1 s k -u∈[0,t) E Y s-u 1 Y t-s ∈A\∂ 1 t-s∈[0,t)

0 h(Yv) dv | H s-u = E a 1 Y s +1 -u ∈A +1 ,...,Y s k -u ∈A k e -s-u 0 h(Yv) dv E Y s-u 1 Y t-s ∈A\∂ e -t-s 0 h(Yv) dv = E a 1

  Xs 1 ∈A 1 ,...,Xs ∈A , s i ∈[τ n i -1 ,τn i ),τm≤s × f (A m , τ m ) ,and then, by definition of the law of (A 0 , Θ 1 , • • • , A m-1 , Θ m , A m , . . .), we deduce thatP µ 0 (X s 1 ∈ A 1 , . . . , X s k ∈ A k , s i ∈ [τ n i -1 , τ n i ), X s+t ∈ A, s + t ∈ [τ n , τ n+1 )) = E µ 0 1 Xs 1 ∈A 1 ,...,Xs k ∈A k , s

i ∈[τ n i -1 ,τn i ) P Xs (X t-s ∈ A, t -s ∈ [0, τ 1 )) . (62)

  1 , • • • , A m-1 , Θ m , A m , . . .), we deduce that P µ 0 (X s 1 ∈ A 1 , . . . , X s k ∈ A k , s i ∈ [τ n i -1 , τ n i ), X s+t ∈ A, s + t ∈ [τ n , τ n+1 )) = E µ 0 1 Xs 1 ∈A 1 ,...,Xs k ∈A k , s i ∈[τ n i -1 ,τn i ) P Xs (X t-s ∈ A, t -s ∈ [τ n-m , τ n-m+1 )) . (63)

According to Ionescu Tulcea's extension Theorem, for any initial distribution µ 0 on F 0 := ∪ s∈P f (N) {s} × (E ∪ ∂) s , there exists a probability measure P µ 0 on Ω = F 0 × F N , such that, for all n ≥ 2,

Then, we can define, for any ω = (a 0 , θ 1 , . . .), the sequence of times (τ n ) n∈N as

and, for all i ∈ St (ω),

, defined on the probability space ( Ω, P) and taking its values in F 0 , is then the size constrained branching process, which concludes the construction of the process.

Remark 26. Observe that the particle system is well defined without Assumptions 1 and 2. However, if Assumption 1 holds true, then ( St , (X i t ) i∈ St ) t≥0 takes its values in

{s} ⊗ E s almost surely. Also, if Assumption 2 holds true, then St is non-empty almost surely, for all t ∈ [0, +∞).

The above construction makes use of the fact that X is a measurable adapted process. Using the fact that it is also assumed progressively measurable, we are in position to prove the following.

Proposition 27. The process ( St , (X i t ) i∈ St ) t≥0 is a Markov process with respect to its natural filtration.

Proof. Fix t ≥ s ≥ 0, and A a measurable subset of F 0 . For all k ≥ 1, all s 1 < s 2 < • • • < s k = s and all measurable subsets A 1 , . . . , A k of F 0 , we have, setting X t := ( St , (X i t ) i∈ St ),

Then, conditioning with respect to G m := σ(Θ 1 , . . . , Θ m , A m ), we have

where f (a, u) := P a X s +1 -u ∈ A +1 , . . . , X s k -u ∈ A k , s i -u ∈ [0, τ 1 ), X t-u ∈ A, t -u ∈ [τ n-m , τ n-m+1 ) .

There are now two possible situations, either m = n or m < n, which we consider separately.