# Binary branching processes with Moran type interactions

A. M. G. Cox111Department of Mathematical Sciences, University of Bath,
Claverton Down, Bath, BA2 7AY, UK. Email<EMAIL_ADDRESS>222Supported by
EPSRC Grant EP/P009220/1, E. Horton333Inria Bordeaux Research Center, 200
Avenue de la Vieille Tour, 33405 Talence, France. Email<EMAIL_ADDRESS>D. Villemonais444Institut Élie Cartan de Lorraine, Bureau 123, Université de
Lorraine, 54506, Vandoeuvre-lès-Nancy Cedex, France. Email:
<EMAIL_ADDRESS>
###### Abstract

The aim of this paper is to study the large population limit of a binary
branching particle system with Moran type interactions: we introduce a new
model where particles evolve, reproduce and die independently and, with a
probability that may depend on the configuration of the whole system, the
death of a particle may trigger the reproduction of another particle, while a
branching event may trigger the death of an other one. We study the occupation
measure of the new model, explicitly relating it to the Feynman-Kac semigroup
of the underlying Markov evolution and quantifying the $\displaystyle L^{2}$
distance between their normalisations. This model extends the fixed size Moran
type interacting particle system discussed in [18, 19, 6, 7, 57] and we will
indeed show that our model outperforms the latter when used to approximate a
birth and death process. We discuss several other applications of our model
including the neutron transport equation [36, 15] and population size
dynamics.

Keywords : interacting particle systems, branching processes, many-to-one,
Markov processes, birth-and-death process, Moran model.

MSC: 82C22, 82C80, 65C05, 60J25, 92D25, 60J80.

##  1 Introduction

Branching processes naturally arise as pertinent models in population size
dynamics [40, 41, 43], neutron transport [14], genetic dynamics [47], growth-
fragmentation processes [4, 3] and cell proliferation kinetics [60], and as
theoretical objects in their own right [40, 17, 37, 38, 39]. 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 [48], the Moran model describes the evolution of
$\displaystyle 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 [28] 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
[12, 21, 20, 57] 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 $\displaystyle
N$ particles, each particle evolves according to an independent copy of a
given Markov process, $\displaystyle 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 $\displaystyle X$. Indeed, letting $\displaystyle
m_{T}$ denote the empirical distribution of the particle system at time
$\displaystyle T$, we show that for any $\displaystyle T\geq 0$,

$\mathbb{E}\left[\Pi_{T}^{A}\,\Pi_{T}^{B}m_{T}f\right]={\mathbf{E}}_{m_{0}}\left[f(X_{t})\exp\left(\int_{0}^{t}b(X_{s})-\kappa(X_{s})\,{\rm
d}s\right)\mathbf{1}_{t<\tau_{\partial}}\right],$ (1)

where $\displaystyle b$ is the branching rate, $\displaystyle\kappa$ is the
soft killing rate, $\displaystyle\tau_{\partial}$ is the absorption time,
$\displaystyle\Pi_{T}^{A}$ and $\displaystyle\Pi_{T}^{B}$ are weights that
compensate for the resampling and selection events that occur up to time
$\displaystyle T$, which we state explicitly in Theorem 6, and
$\displaystyle{\mathbf{E}}_{\cdot}$ is the expectation with respect to the law
of $\displaystyle X$. As we will later discuss, by choosing particular model
parameters, we recover both the classical many-to-one formula [35] and the
unbiased estimator proved in [57] for fixed size Moran type genetic
algorithms. The second part of Theorem 6 then provides a precise bound on the
$\displaystyle 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 $\displaystyle\\{N_{min},\ldots,N_{max}\\}$,
where $\displaystyle 2\leq N_{min}\leq N_{max}<+\infty$. In order to constrain
the size of the process, when the size of the population reaches
$\displaystyle N_{max}$ (resp. $\displaystyle N_{min}$) and a natural
branching (resp. killing) event occurs, we set the probability of selection
(resp. resampling) to be $\displaystyle 1$. We often refer to this model as
the $\displaystyle N_{min}$–$\displaystyle N_{max}$ model. As will be clear in
the next section, the choice $\displaystyle N_{min}=0$ and $\displaystyle
N_{max}=+\infty$ 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 $\displaystyle N_{min}=N_{max}$) to branching
population models (when $\displaystyle N_{min}=0$ and $\displaystyle
N_{max}=+\infty$) with or without a selection/resampling mechanism.

For comparable results concerning the resampling algorithm, we refer the
reader to [29, 51, 13, 49] and to [27, 22, 45] 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 [18, 19] 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 [6, 7, 32, 51, 56, 12] 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 [52], 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 [59] to allow for random control functions. We refer the reader to
[52, 31, 55, 58] for discrete time versions of this process and to [30] 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, [58]. We refer the
reader to [54, 42, 50, 46, 44, 2, 1] 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 $\displaystyle N_{min}$–$\displaystyle 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 $\displaystyle
N_{min}$–$\displaystyle 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 $\displaystyle 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 $\displaystyle 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
$\displaystyle\sqrt{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 [36, 15], we
overcome this difficulty using the notion of $\displaystyle h$-transform
introduced in [15] for the NTE. We show that particles whose dynamics are
given by an appropriately $\displaystyle h$-transformed process fulfil the
assumptions of our main theorem, and thus, we may use the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ model with the $\displaystyle 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
$\displaystyle N_{min}$–$\displaystyle N_{max}$ model and the fixed size Moran
type interacting particle system (IPS) for the approximation of quantities
related to non-conservative semigroups (cf. [18, 19, 6, 7]), which corresponds
to the particular case $\displaystyle N_{min}=N_{max}$ with no branching. We
observe that in the case of a birth-and-death process, $\displaystyle X$, with
branching rate given by $\displaystyle X_{t}\wedge M$ for some $\displaystyle
M>0$ and state-dependent Markov transitions, the computational cost, measured
as the number of resampling/selection events, is significantly lower when
using the $\displaystyle N_{min}$–$\displaystyle N_{max}$ algorithm. We also
study the bias and standard deviation of an estimator for the normalised
empirical stationary distribution for each of the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ and fixed size Moran type IPS. Letting
$\displaystyle M$ go to infinity, we provide numerical evidence that the
number of interactions, bias and standard deviation stabilise for the
$\displaystyle N_{min}$–$\displaystyle 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.

##  2 Main results

###  2.1 Description of the model

Let $\displaystyle(\Omega,\mathcal{F},(X_{t})_{t\in[0,+\infty)})$ be a
continuous time progressively measurable Markov process with values in a
measurable state space $\displaystyle E\cup{\partial}$, where
$\displaystyle{\partial}$ is an absorbing (measurable) set such that
$\displaystyle{\partial}\cap E=\emptyset$. Denoting by
$\displaystyle\tau_{\partial}=\inf\\{t\geq 0,\ X_{t}\in{\partial}\\}$ its
absorption time, it follows that

$\displaystyle X_{t}\in{\partial},\ \forall t\geq\tau_{\partial}.$

We denote by $\displaystyle{\mathbf{P}}_{x}$ its law when initiated at
$\displaystyle x\in E$ and by $\displaystyle{\mathbf{E}}_{x}$ the
corresponding expectation operator. We extend, whenever necessary, any
measurable function $\displaystyle f:E\to[0,+\infty)$ by $\displaystyle
f\equiv 0$ on $\displaystyle{\partial}$. We also assume that we are given two
bounded functions $\displaystyle b:E\to\mathbb{R}_{+}$ and
$\displaystyle\kappa:E\to\mathbb{R}_{+}$, denoting.

We will shortly define a particle system, where particles move according to
copies of $\displaystyle 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 $\displaystyle X$ are absorbed simultaneously
with probability $\displaystyle 0$.

Assumption 1. For any $\displaystyle x\in E$ and $\displaystyle
t\in[0,+\infty)$, $\displaystyle{\mathbf{P}}_{x}(\tau_{\partial}=t)=0$ and
$\displaystyle{\mathbf{P}}_{x}(\tau_{\partial}>t)>0$.

We also introduce the following notation. Denote by
$\displaystyle\mathcal{P}_{f}(\mathbb{N})$ the collection of finite subsets of
$\displaystyle\mathbb{N}:=\\{1,2,\dots\\}$. Then, for $\displaystyle
i\in\mathbb{N}$, we let

$\displaystyle\displaystyle
b_{i}:(E\cup{\partial})^{\mathcal{P}_{f}(\mathbb{N})}\to[0,+\infty),$
$\displaystyle\displaystyle\kappa_{i}:(E\cup{\partial})^{\mathcal{P}_{f}(\mathbb{N})}\to[0,+\infty)$

be a collection of bounded measurable functions. Also for $\displaystyle
i\in\mathbb{N}$, let

$\displaystyle\displaystyle
p_{i}:(E\cup{\partial})^{\mathcal{P}_{f}(\mathbb{N})}\to[0,1],$
$\displaystyle\displaystyle
q_{i}:(E\cup{\partial})^{\mathcal{P}_{f}(\mathbb{N})}\to[0,1]$

be measurable functions. We will assume that, for any $\displaystyle
s\in\mathcal{P}_{f}(\mathbb{N})$, for any collection of points
$\displaystyle\\{x_{i}:i\in s,x_{i}\in E\\}$, and any $\displaystyle i_{0}\in
s$,

$\displaystyle\displaystyle b^{i_{0}}(x_{i},i\in s)-\kappa^{i_{0}}(x_{i},i\in
s)=b(x_{i_{0}})-\kappa(x_{i_{0}})$ (2)

and that

$\displaystyle\displaystyle p^{i_{0}}(x_{i},i\in s)=0\text{ whenever }|s|=1,$
(3)

where $\displaystyle|s|$ denotes the number of elements in the set
$\displaystyle 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 $\displaystyle\bar{N}_{0}\geq 2$. We
consider the particle system
$\displaystyle(\bar{S}_{t},(X^{i}_{t})_{i\in\bar{S}_{t}})_{t\in[0,+\infty)}$,
where the component $\displaystyle\bar{S}_{t}\in\mathcal{P}_{f}(\mathbb{N})$
is the set enumerating the particles in the system at time $\displaystyle t$,
initiated with $\displaystyle\bar{S}_{0}=\\{1,2,\ldots,\bar{N}_{0}\\}$ and
$\displaystyle X^{i}_{0}=x_{i}\in E$ for all $\displaystyle i\in\bar{S}_{0}$.
The number of particles in the system at any time $\displaystyle t$ is denoted
by $\displaystyle\bar{N}_{t}=|\bar{S}_{t}|$.

Evolution of the BBMMI.

  1. 1.

The particles $\displaystyle X^{i}$, $\displaystyle i\in\bar{S}_{0}$, evolve
as independent copies of $\displaystyle X$, and we consider the following
times:

$\displaystyle\displaystyle\tau^{b,i}_{1}:=\inf\\{t\geq 0,\
\int_{0}^{t}b^{i}(X^{i}_{s},\,i\in\bar{S}_{0})\,ds\geq e_{1}^{i}\\},$

and

$\displaystyle\displaystyle\tau^{\kappa,i}_{1}:=\inf\\{t\geq 0,\
\int_{0}^{t}\kappa^{i}(X^{i}_{s},\,i\in\bar{S}_{0})\,ds\geq E_{1}^{i}\\}$

and

$\displaystyle\displaystyle\tau^{{\partial},i}_{1}:=\inf\\{t\geq 0,\
X^{i}_{t}\in{\partial}\\},$

where $\displaystyle e_{1}^{i},E_{1}^{i}$, $\displaystyle
i=1,\dots,\bar{N}_{0}$ are exponential random variables with parameter
$\displaystyle 1$, and are independent of each other and everything else.

  2. 2.

Denoting by $\displaystyle i_{0}$ the index of the (unique) particle such that
$\displaystyle\tau^{b,i_{0}}_{1}\wedge\tau^{\kappa,i_{0}}_{1}\wedge\tau^{\partial,i_{0}}_{1}=\tau_{1}$,
where
$\displaystyle\tau_{1}=\min_{i\in\bar{S}_{0}}\tau^{b,i}_{1}\wedge\tau^{\kappa,i}_{1}\wedge\tau^{\partial,i}_{1}$,
we set $\displaystyle\bar{S}_{t}=\bar{S}_{0}$ for all $\displaystyle
t<\tau_{1}$ and

    1. (a)

if $\displaystyle\tau_{1}=\tau^{b,i_{0}}_{1}$, then a branching event occurs;

    2. (b)

if $\displaystyle\tau_{1}=\tau^{\kappa,i}_{1}$, then a soft killing event
occurs;

    3. (c)

if $\displaystyle\tau_{1}=\tau^{{\partial},i_{0}}_{1}$, then a hard killing
event occurs.

  3. 3.

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

    1. killing: if a (hard or soft) killing event occurred at the preceding step, then we say that $\displaystyle i_{0}$ is killed at time $\displaystyle\tau_{1}$ and

       * •

with probability $\displaystyle
p^{i_{0}}(X^{i}_{\tau_{1}},\,i\in\bar{S}_{0})$, the particle $\displaystyle
i_{0}$ is removed from the system and a resampling event occurs: one chooses
$\displaystyle j_{0}$ uniformly in $\displaystyle S_{0}\setminus\\{i_{0}\\}$
and sets

$\displaystyle X^{\max S_{0}+1}_{\tau_{1}}=X^{j_{0}}_{\tau_{1}}\text{ and
}\bar{S}_{\tau_{1}}:=\bar{S}_{0}\setminus\\{i_{0}\\}\cup\\{\max\bar{S}_{0}+1\\};$

observe that the number of particles in the system at time
$\displaystyle\tau_{1}$ is then $\displaystyle\bar{N}_{\tau_{1}}=\bar{N}_{0}$;
we say that $\displaystyle j_{0}$ is duplicated at time
$\displaystyle\tau_{1}$.

       * •

with probability $\displaystyle
1-p^{i_{0}}(X^{i}_{\tau_{1}},\,i\in\bar{S}_{0})$, the particle $\displaystyle
i_{0}$ is removed from the system; the set of particles at time
$\displaystyle\tau_{1}$ is enumerated by
$\displaystyle\bar{S}_{\tau_{1}}:=\bar{S}_{0}\setminus\\{i_{0}\\}$ and the
number of particles in the system is then
$\displaystyle\bar{N}_{\tau_{1}}=\bar{N}_{0}-1$;

    2. branching: if a branching event occurred at the preceding step, then we say that $\displaystyle i_{0}$ has branched at time $\displaystyle\tau_{1}$ and

       * •

with probability $\displaystyle
q^{i_{0}}(X^{i}_{\tau_{1}},\,i\in\bar{S}_{0})$, a new particle is added to the
system at position $\displaystyle X^{i_{0}}_{\tau_{1}}$ and a selection event
occurs: one chooses $\displaystyle j_{0}$ at random uniformly in
$\displaystyle S_{0}\cup\\{\max S_{0}+1\\}$ and removes the particle
$\displaystyle j_{0}$ from the system:

$\displaystyle X^{\max S_{0}+1}_{\tau_{1}}=X^{i_{0}}_{\tau_{1}}\text{ and
}\bar{S}_{\tau_{1}}:=\\{\max\bar{S}_{0}+1\\}\cup\bar{S}_{0}\setminus\\{j_{0}\\};$

observe that the number of particles at time $\displaystyle\tau_{1}$ is thus
$\displaystyle\bar{N}_{\tau_{1}}=\bar{N}_{0}$; we say that $\displaystyle
j_{0}$ is removed at time $\displaystyle\tau_{1}$ and that $\displaystyle\max
S_{0}+1$ is born at time $\displaystyle\tau_{1}$;

       * •

with probability $\displaystyle
1-q^{i_{0}}(X^{i}_{\tau_{1}},\,i\in\bar{S}_{0})$, a new particle is added to
the system at position $\displaystyle X^{i_{0}}_{\tau_{1}}$:

$\displaystyle X^{\max S_{0}+1}_{\tau_{1}}=X^{i_{0}}_{\tau_{1}}\text{ and
}\bar{S}_{\tau_{1}}=\bar{S}_{0}\cup\\{\max\bar{S}_{0}\\},$

and we say that $\displaystyle\max S_{0}+1$ is born at time
$\displaystyle\tau_{1}$.

After time $\displaystyle\tau_{1}$ the system evolves as independent copies of
$\displaystyle X$ until the next killing/branching event, denoted by
$\displaystyle\tau_{2}$, and at time $\displaystyle\tau_{2}$ it may undergo a
resampling/selection event as above. By iteration, we define the sequence
$\displaystyle\tau_{0}:=0<\tau_{1}<\tau_{2}<\cdots<\tau_{n}<\cdots$.

(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.

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

We will also make use of the following assumption, which ensures that the
process described above is well defined at any time $\displaystyle t\geq 0$.

Assumption 2. The sequence $\displaystyle(\tau_{n})_{n\in\mathbb{N}}$
converges to $\displaystyle+\infty$ 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 $\displaystyle b^{i_{0}}$ and $\displaystyle\kappa^{i_{0}}$ are
the rate at which the particle $\displaystyle i$ in the system branches and is
killed respectively (in addition to hard killing events occurring when a
particle hits the absorbing set $\displaystyle\partial$). The simplest
situation is when $\displaystyle b^{i_{0}}(x^{i},\,i\in s)=b(x_{i_{0}})$ and
$\displaystyle\kappa^{i_{0}}(x_{i},\,i\in s)=\kappa(x_{i_{0}})$, but other
configurations are possible. We emphasise that the balance condition (2) is
crucial to ensure the relation with the semigroup $\displaystyle Q$ defined in
the next section, and should not be considered a technical assumption.

######  Remark 3.

The functions $\displaystyle p^{i_{0}}$ (respectively $\displaystyle
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 $\displaystyle 0$ particles. We emphasize that, since
$\displaystyle\partial$ is a set, $\displaystyle 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 $\displaystyle N_{min}\neq 1$ and $\displaystyle N_{max}$,
with $\displaystyle 0\leq N_{min}\leq N_{max}$. In order to do so, one simply
chooses

$\displaystyle p^{i_{0}}(x_{i},\,i\in s)=\mathbf{1}_{\\{|s|=N_{min}\\}}\text{
and }q^{i_{0}}(x_{i},\,i\in s)=\mathbf{1}_{\\{|s|=N_{max}\\}}.$

In what follows, we refer to this dynamic as the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ process. Other natural choices include, for
instance,

$\displaystyle p^{i_{0}}(x_{i},\,i\in s)=\frac{1}{\sharp
s+1}\mathbf{1}_{|s|\geq 2}\text{ and }q^{i_{0}}(x_{i},\,i\in
s)=1-\frac{1}{\sharp 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 $\displaystyle X$ is only assumed to be
Markov, it can carry the time as a component, which would allow $\displaystyle
N_{min}$ and $\displaystyle 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 $\displaystyle p,q$ depend on time.

######  Remark 5.

In our model, we consider general killing, which corresponds to a conservative
Markov process $\displaystyle 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
$\displaystyle\kappa$ 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 $\displaystyle 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 [7, 32, 56, 57] to our
setting.

In the rest of this article, we denote by $\displaystyle\sigma_{n}$ and
$\displaystyle\rho_{n}$ the $\displaystyle n^{th}$ selection and resampling
events, respectively. We let $\displaystyle A_{t}$ denote the total number of
resampling events up to time $\displaystyle t$, $\displaystyle B_{t}$ denote
the total number of selection events up to time $\displaystyle t$ and we set
$\displaystyle C_{t}=\inf\\{n\geq 0:t<\tau_{n+1}\\}$ to be the total number of
events up to time $\displaystyle t$. We will also use the following quantities
throughout,

$\Pi_{t}^{A}\coloneqq\prod_{n=1}^{A_{t}}\left(\frac{\bar{N}_{\rho_{n}}-1}{\bar{N}_{\rho_{n}}}\right),\qquad\Pi_{t}^{B}\coloneqq\prod_{i=1}^{B_{t}}\left(\frac{\bar{N}_{\sigma_{n}}+1}{\bar{N}_{\sigma_{n}}}\right).$
(4)

###  2.2 $\displaystyle L^{2}$ bounds for the empirical distribution of the
process

For all bounded measurable functions $\displaystyle f:E\to\mathbb{R}$, all
$\displaystyle t\geq 0$ and all $\displaystyle x\in E$, we set

$\displaystyle\displaystyle Q_{t}f(x)$
$\displaystyle\displaystyle={\mathbf{E}}_{x}\left[f(X_{t})\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,{\rm
d}s\right)\mathbf{1}_{t<\tau_{\partial}}\right],$

where we recall that $\displaystyle\tau_{\partial}=\inf\\{t\geq
0:X_{t}\in\partial\\}$. This defines a Feynman-Kac semigroup
$\displaystyle(Q_{t})_{t\geq 0}$, which is related to the binary branching
model where particles move as copies of $\displaystyle X$ that are killed at
rate $\displaystyle\kappa$ and branch at rate $\displaystyle b$ resulting in
the creation of two independent copies of the original particle. The relation
between $\displaystyle Q$ and this process is given by the well known many-to-
one formula (see for instance [35] and references therein, see also Remark 8).
In this section, we consider the BBMMI defined above and study its relation
with $\displaystyle Q$.

In what follows, we let $\displaystyle m_{t}$ (resp
$\displaystyle\hat{m}_{t}$) denote the occupation measure (resp. normalised
occupation measure) of the particle system,

$m_{t}=\sum_{i\in\bar{S}_{t}}\delta_{X_{t}^{i}},\qquad t\geq
0,\quad\hat{m}_{t}=\frac{1}{\bar{N}_{0}}m_{t}.$ (5)

The following result first provides an equality between a penalised version of
the occupation measure of the system at any time $\displaystyle T$ and the
semigroup $\displaystyle Q$ at time $\displaystyle T$. The second part is
dedicated to the control of the $\displaystyle L^{2}$ distance between
renormalised versions of the semigroup $\displaystyle Q$ and the occupation
measure of the process.

######  Theorem 6.

Under Assumptions 1 and 2, the BBMMI satisfies, for all time $\displaystyle
T\geq 0$ and all bounded measurable functions $\displaystyle
f:E\to\mathbb{R}$, the following many-to-one formula

$m_{0}Q_{T}f=\mathbb{E}\left(\Pi^{A}_{T}\,\Pi^{B}_{T}m_{T}f\right).$ (6)

Moreover,

$\left\|\frac{m_{0}Q_{T}f}{m_{0}Q_{T}\mathbf{1}_{E}}-\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}\mathbf{1}_{m_{T}\neq
0}\right\|_{2}\leq
C\,\exp(c\|b\|_{\infty}T)\frac{\|f\|_{\infty}}{m_{0}Q_{T}\mathbf{1}_{E}/{\bar{N}_{0}}}\,\frac{1}{\sqrt{\bar{N}_{0}}},$
(7)

where $\displaystyle 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, $\displaystyle m_{0}$ is assumed deterministic with total mass
$\displaystyle\bar{N}_{0}$. The result immediately extends to the case where
$\displaystyle 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,
$\displaystyle\hat{m}$, in which case (6) becomes
$\displaystyle\hat{m}_{0}Q_{T}f=\mathbb{E}\left(\Pi^{A}_{T}\,\Pi^{B}_{T}\hat{m}_{T}f\right)$
and (7) can be written as

$\left\|\frac{\hat{m}_{0}Q_{T}f}{\hat{m}_{0}Q_{T}\mathbf{1}_{E}}-\frac{\hat{m}_{T}(f)}{\hat{m}_{T}(\mathbf{1}_{E})}\mathbf{1}_{\hat{m}_{T}\neq
0}\right\|_{2}\leq
C\,\exp(c\|b\|_{\infty}T)\frac{\|f\|_{\infty}}{\hat{m}_{0}Q_{T}\mathbf{1}_{E}}\,\frac{1}{\sqrt{\bar{N}_{0}}},$

making explicit the ‘classical’ dependence $\displaystyle{\bar{N}_{0}}^{-1/2}$
on the number of particles which appears on the right-hand side.

######  Remark 8.

When $\displaystyle p=q\equiv 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 $\displaystyle\Pi^{A}_{T}=\Pi^{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 $\displaystyle p=q\equiv 1$ (so that the size of the process remains
constant equal to $\displaystyle\bar{N}_{0}$), Theorem 6 extends the results
of [57] by allowing selection events. Similar results were also obtained in
[51] (see also references therein and the recent paper [12] for an adaptation
to discrete state spaces, with more general interaction mechanisms), where the
process is constrained to remain of constant size $\displaystyle\bar{N}_{0}$
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 $\displaystyle Q$ (see
the above references and also [26, 49, 53] 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
$\displaystyle(\xi_{t})_{t\geq 0}$ from their point of creation, and when at
position $\displaystyle x\in E$, at rate $\displaystyle\beta(x)$, the particle
is removed from the system and replaced by two particles at (random) positions
$\displaystyle x_{1},x_{2}\in E$, which may be different to $\displaystyle x$.
Particles also undergo (soft) killing at rate $\displaystyle\kappa(x)$.
Denoting by $\displaystyle(\psi_{t})_{t\geq 0}$ the linear semigroup of this
process, in this case, the many-to-one formula is given by

$\displaystyle\psi_{t}[g](x)=\mathbf{E}_{x}\left[{\rm
e}^{\int_{0}^{t}b(X_{s})-\kappa(X_{s})\mathrm{d}s}g(X_{t})\mathbf{1}_{t<\tau_{\partial}}\right],$

where $\displaystyle(X_{t})_{t\geq 0}$ behaves like
$\displaystyle(\xi_{t})_{t\geq 0}$, except that at rate $\displaystyle
2b(X_{t})$, the process $\displaystyle X_{t}$ jumps to a new location given by
$\displaystyle(\delta_{x_{1}}+\delta_{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
$\displaystyle(\bar{S}_{t},(X^{i}_{t})_{t\in\bar{S}_{t}})_{t\geq 0}$ that is
used in the proof of Theorem 6 is the identity

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}\left[Q_{T-(t+s)}f(X^{1}_{t+s})\mathbf{1}_{t+s<\tau_{1}}\right]$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[Q_{T-(t+s)}f(Y^{1}_{t+s})e^{-\int_{t}^{t+s}h_{u}\,{\rm
d}u}\mathbf{1}_{t+s<\tau_{\partial}}\right]\mathbf{1}_{t<\tau_{1}},$ (8)

where $\displaystyle\texttt{E}^{t}$ is the expectation on a probability space
$\displaystyle\Omega^{\prime}$ such that $\displaystyle(Y^{i}_{u})_{u\geq t}$,
$\displaystyle i\in\bar{S}_{0}$, are independent copies of $\displaystyle X$,
starting from $\displaystyle X^{i}_{t}$ at time $\displaystyle t$, and where

$\displaystyle
h_{u}:=\sum_{j\in\bar{S}_{0}}b^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}})+\kappa^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}}).$

This identity is derived from the Markov property, and the assumption on the
dynamics of the process
$\displaystyle(\bar{S}_{t},(X^{i}_{t})_{t\in\bar{S}_{t}})_{t\geq 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 $\displaystyle(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq
0},\mathbb{P})$ and a progressively measurable process
$\displaystyle(\bar{S}_{t},(X^{i}_{t})_{t\in\bar{S}_{t}})_{t\geq 0}$ defined
on this space, such that the rates $\displaystyle
b^{i},\kappa^{i},p^{i},q^{i}$ are no longer functions, but merely
progressively measurable processes, satisfying the equivalents555Technically,
we need that the integrals of the respective processes are equal up to
modification, that is, if we define $\displaystyle
I^{i,1}_{t}:=\int_{0}^{t}(b_{s}^{i}-\kappa_{s}^{i})\,{\rm d}s$ and
$\displaystyle
I^{i,2}_{t}:=\int_{0}^{t}(b^{i}(X_{s}^{i})-\kappa(X_{s}^{i}))\,{\rm d}s$, then
we require the processes $\displaystyle I^{i,1}$ and $\displaystyle I^{i,2}$
to agree up to a modification. Note that as a consequence of the difference
$\displaystyle b(\cdot)-\kappa(\cdot)$ being bounded, this implies that the
processes are also indistinguishable. of (2), (3):

$\displaystyle\displaystyle b^{i}_{t}-\kappa^{i}_{t}$
$\displaystyle\displaystyle=b(X_{t}^{i})-\kappa(X_{t}^{i})$
$\displaystyle\displaystyle p^{i}_{t}$
$\displaystyle\displaystyle=0\quad\text{ whenever }|\bar{S}_{t}|=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
$\displaystyle\mathbb{E}^{t}_{x}$ as the expectation conditional on the full
sigma-algebra $\displaystyle\mathcal{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 $\displaystyle 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 $\displaystyle[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 $\displaystyle b^{i}$ and $\displaystyle\kappa^{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
$\displaystyle(\mathcal{F}_{t})_{t\geq 0}$-progressively measurable, for an
appropriate filtration, which may then be larger than the filtration generated
by the particle histories.

##  3 Applications

###  3.1 Population size dynamics with constrained population size

In this section, we introduce a branching model for a population where the
size is constrained between $\displaystyle N_{min}\geq 2$ and $\displaystyle
N_{max}\geq 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 $\displaystyle x\in E:=\mathbb{Z}_{+}^{d}$, $\displaystyle d\geq
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

$\displaystyle\begin{cases}\beta_{i}(x)&=\text{ transition rate from
$\displaystyle x$ to $\displaystyle x+e_{i}$},\\\ \delta_{i}(x)&=\text{
transition rate from $\displaystyle x$ to $\displaystyle
x-e_{i}$},\end{cases}$

where $\displaystyle e_{i}\in\mathbb{Z}_{+}^{d}$ denotes the vector with
$\displaystyle 1$ in the $\displaystyle i^{th}$ position and $\displaystyle 0$
everywhere else,
$\displaystyle\beta_{i},\delta_{i}:\mathbb{Z}_{+}^{d}\to[0,+\infty)$ for all
$\displaystyle i\in\\{1,\ldots,d\\}$, and $\displaystyle\delta_{i}(x)=0$ for
all $\displaystyle x=(x_{1},\ldots,x_{d})$ such that $\displaystyle x_{i}=0$
(so that the process cannot leave the set $\displaystyle\mathbb{Z}_{+}^{d}$).

When its health state is $\displaystyle x\in\mathbb{Z}_{+}^{d}$, at rate
$\displaystyle\kappa(x)$ an individual dies, where
$\displaystyle\kappa:\mathbb{Z}_{+}^{d}\to[0,+\infty)$, and at rate
$\displaystyle b(x)$, it gives birth to an identical individual with equal
health state $\displaystyle x$, where $\displaystyle
b:\mathbb{Z}_{+}^{d}\to[0,+\infty)$ 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 $\displaystyle 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 $\displaystyle
N_{max}$ by removing an individual chosen randomly and uniformly in the
population when the population size reaches $\displaystyle N_{max}+1$, and the
total number of individuals is lower bounded by $\displaystyle N_{min}$ by
cloning an individual chosen randomly and uniformly when the population size
reaches $\displaystyle N_{min}-1$ (in the situation where $\displaystyle
N_{min}=0$, this last mechanism does not operate).

Formally, the state space of our population model is $\displaystyle
F:=\bigcup_{N=N_{min}}^{N_{max}}[E]^{N}$, where $\displaystyle[E]^{N}$ denotes
the unordered $\displaystyle N$-tuples of $\displaystyle E$. We represent the
elements of $\displaystyle F$ by
$\displaystyle{\mathbf{x}}=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{N}]$ and
note that repetitions are allowed. The extended infinitesimal generator of the
process acting on bounded mesurable functions $\displaystyle
f:F\to\mathbb{R}$, is given, for any
$\displaystyle{\mathbf{x}}=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{N}]\in F$,
by

$\displaystyle\displaystyle\bar{L}f({\mathbf{x}})$
$\displaystyle\displaystyle=\sum_{k=1}^{N}\sum_{i=1}^{d}\left(\beta_{i}({\mathbf{x}}_{k,i})\Delta_{{\mathbf{x}},{\mathbf{x}}+e^{N}_{k,i}}f+\delta_{i}({\mathbf{x}}_{k,i})\Delta_{{\mathbf{x}},{\mathbf{x}}-e^{N}_{k,i}}f\right)$
$\displaystyle\displaystyle\qquad\qquad+\sum_{k=1}^{N}b({\mathbf{x}}_{k})\Delta_{{\mathbf{x}},{\mathbf{x}}\cup[{\mathbf{x}}_{k}]}f+\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})\Delta_{{\mathbf{x}},{\mathbf{x}}\setminus[{\mathbf{x}}_{k}]}f({\mathbf{x}})$
$\displaystyle\displaystyle\qquad\qquad+\mathbf{1}_{N=N_{min}}\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})\frac{1}{N-1}\sum_{\ell=1,\ell\neq
k}^{N}\Delta_{{\mathbf{x}}\setminus[{\mathbf{x}}_{k}],{\mathbf{x}}\cup[{\mathbf{x}}_{\ell}]\setminus[{\mathbf{x}}_{k}]}f$
$\displaystyle\displaystyle\qquad\qquad+\mathbf{1}_{N=N_{max}}\sum_{k=1}^{N}b({\mathbf{x}}_{k})\frac{1}{N+1}\sum_{\ell=1}^{N}(1+\mathbf{1}_{\ell=k})\Delta_{{\mathbf{x}}\cup[{\mathbf{x}}_{k}],{\mathbf{x}}\cup[{\mathbf{x}}_{k}]\setminus[{\mathbf{x}}_{\ell}]}f,$

where

  * •

$\displaystyle{\mathbf{x}}_{k,i}$ is the $\displaystyle i^{\text{th}}$
component of $\displaystyle{\mathbf{x}}_{k}$,

  * •

$\displaystyle{\mathbf{x}}+e^{N}_{k,i}$ is $\displaystyle{\mathbf{x}}$ but
with $\displaystyle{\mathbf{x}}_{k,i}$ replaced by
$\displaystyle{\mathbf{x}}_{k,i}+1$,

  * •

$\displaystyle{\mathbf{x}}-e^{N}_{k,i}$ is $\displaystyle{\mathbf{x}}$ but
with $\displaystyle{\mathbf{x}}_{k,i}$ replaced by
$\displaystyle{\mathbf{x}}_{k,i}-1$,

  * •

for all $\displaystyle\mathbf{x},\mathbf{y}$,
$\displaystyle\Delta_{{\mathbf{x}},\mathbf{y}}f:=f(\mathbf{y})\mathbf{1}_{\mathbf{y}\in
F}-f({\mathbf{x}})\mathbf{1}_{\mathbf{x}\in F}$,

  * •

for all $\displaystyle y\in E$,
$\displaystyle{\mathbf{x}}\cup[y]=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{N},y]$,

  * •

$\displaystyle{\mathbf{x}}\setminus[{\mathbf{x}}_{k}]=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{k-1},{\mathbf{x}}_{k+1},\ldots,{\mathbf{x}}_{N}]$.

In a similar manner to the previous section, we denote by $\displaystyle
A_{T}$ the number of times an individual is duplicated by the size
constraining mechanism before time $\displaystyle T$, and by $\displaystyle
B_{T}$ the number of times an individual is removed from the population by the
size constraining mechanism. We also set
$\displaystyle({\mathbf{X}}(t))_{t\geq 0}$ to be the population process,
$\displaystyle(|{\mathbf{X}}(t)|)_{t\geq 0}$ its size, and
$\displaystyle(X_{t})_{t\in\geq 0}$ a multi-dimensional birth and death
process with (state dependent) transition rates $\displaystyle\beta$ and
$\displaystyle\delta$. The notation $\displaystyle\mathbb{P}_{\mathbf{x}}$
(respectively $\displaystyle\mathbf{P}_{x}$) is used for the law of the
population process $\displaystyle{\mathbf{X}}$ starting from
$\displaystyle{\mathbf{x}}\in F$ (respectively for the law of $\displaystyle
X$ starting from $\displaystyle x\in E$). We denote as usual by
$\displaystyle\mathbb{E}_{x}$ and $\displaystyle\mathbf{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 $\displaystyle\kappa$ 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 $\displaystyle V$ and $\displaystyle\kappa$ are co-
monotone if $\displaystyle(\kappa(x)-\kappa(y))(V(x)-V(y))\geq 0$ for all
$\displaystyle x,y\in E$. We say that $\displaystyle V$ and
$\displaystyle\kappa$ are almost co-monotone if there exists a function
$\displaystyle\kappa^{\prime}$ co-monotone with $\displaystyle V$ such that
$\displaystyle\kappa-\kappa^{\prime}$ is bounded. We say that $\displaystyle
V$ goes to infinity at infinity if the set of $\displaystyle x\in E$ such that
$\displaystyle V(x)<c$ is finite for all $\displaystyle c>0$.

######  Proposition 13.

Assume that there exists a function $\displaystyle V:E\to[1,+\infty)$ almost
co-monotone with $\displaystyle\kappa$, such that $\displaystyle V$ tends to
infinity at infinity, and such that, for all $\displaystyle x\in E$,

$\displaystyle\displaystyle\sum_{i=1}^{d}\left[\beta_{i}(x)(V(x+e_{i})-V(x))+\delta_{i}(x)(V(x-e_{i})-V(x))\right]\leq
CV(x),$ (9)

for some constant $\displaystyle C\in\mathbb{R}$. Then the process with
transition rates given by $\displaystyle\bar{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 $\displaystyle\bar{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 $\displaystyle N_{max}=+\infty$, or if we
assume that, immediately after each branching time and with probability
$\displaystyle 1+\frac{1}{N_{max}+1}-\frac{1}{N}$, a uniformly chosen
individual is removed from the population, and/or that, immediately after each
killing time and with probability $\displaystyle\frac{N_{min}-1}{N}$, a
uniformly chosen individual is duplicated (where $\displaystyle 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

$\displaystyle\displaystyle\bar{V}:$ $\displaystyle\displaystyle
F\to\mathbb{[}1,+\infty)$ $\displaystyle\displaystyle:$
$\displaystyle\displaystyle{\mathbf{x}}=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{N}]\mapsto\sum_{k=1}^{N}V({\mathbf{x}}_{k}),$

so that $\displaystyle\bar{V}$ goes to infinity at infinity and for
$\displaystyle{\mathbf{x}}=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{N}]\in F$,

$\displaystyle\displaystyle\bar{L}\bar{V}({\mathbf{x}})$
$\displaystyle\displaystyle=\sum_{k=1}^{N}\sum_{i=1}^{d}\left[\beta_{i}({\mathbf{x}}_{k})(V({\mathbf{x}}_{k}+e_{i})-V({\mathbf{x}}_{k}))+\delta_{i}({\mathbf{x}}_{k})(V({\mathbf{x}}_{k}-e_{i})-V({\mathbf{x}}_{k}))\right]$
$\displaystyle\displaystyle\qquad\qquad+\sum_{k=1}^{N}\left(b({\mathbf{x}}_{k})-\kappa({\mathbf{x}}_{k})\right)V({\mathbf{x}}_{k})$
$\displaystyle\displaystyle\qquad\qquad+\mathbf{1}_{N=N_{min}}\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})\frac{1}{N-1}\sum_{\ell=1,\ell\neq
k}^{N}V({\mathbf{x}}_{l})$
$\displaystyle\displaystyle\qquad\qquad-\mathbf{1}_{N=N_{max}}\sum_{k=1}^{N}b({\mathbf{x}}_{k})\frac{1}{N+1}\sum_{\ell=1}^{N}(1+\mathbf{1}_{\ell=k})V({\mathbf{x}}_{\ell}).$

Now let $\displaystyle\kappa^{\prime}$ be co-monotone with $\displaystyle V$
such that $\displaystyle\kappa^{\prime\prime}=\kappa-\kappa^{\prime}$ is
bounded. Then, using the fact that $\displaystyle b$ is bounded, we deduce
that

$\displaystyle\displaystyle\bar{L}\bar{V}({\mathbf{x}})$
$\displaystyle\displaystyle\leq(C+\|b\|_{\infty})\overline{V}({\mathbf{x}})+\mathbf{1}_{N=N_{min}}\left(\frac{1}{N-1}\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})\sum_{{\stackrel{{\scriptstyle\ell=1}}{{\ell\neq
k}}}}^{N}V({\mathbf{x}}_{\ell})-\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})V({\mathbf{x}}_{k})\right)$
$\displaystyle\displaystyle=(C+\|b\|_{\infty})\overline{V}({\mathbf{x}})+\frac{\mathbf{1}_{N=N_{min}}}{N-1}\left(\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})\sum_{\ell=1}^{N}V({\mathbf{x}}_{\ell})-N\sum_{k=1}^{N}\kappa({\mathbf{x}}_{k})V({\mathbf{x}}_{k})\right)$
$\displaystyle\displaystyle=(C+\|b\|_{\infty})\overline{V}({\mathbf{x}})+\frac{N^{2}\mathbf{1}_{N=N_{min}}}{N-1}\left[\left(\sum_{k=1}^{N}\frac{\kappa({\mathbf{x}}_{k})}{N}\right)\left(\sum_{k=1}^{N}\frac{V({\mathbf{x}}_{k})}{N}\right)-\sum_{k=1}^{N}\frac{\kappa({\mathbf{x}}_{k})V({\mathbf{x}}_{k})}{N}\right]$
$\displaystyle\displaystyle\leq(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})\overline{V}({\mathbf{x}})$
$\displaystyle\displaystyle\qquad\qquad+\frac{N^{2}\mathbf{1}_{N=N_{min}}}{N-1}\left[\left(\sum_{k=1}^{N}\frac{\kappa^{\prime}({\mathbf{x}}_{k})}{N}\right)\left(\sum_{k=1}^{N}\frac{V({\mathbf{x}}_{k})}{N}\right)-\sum_{k=1}^{N}\frac{\kappa^{\prime}({\mathbf{x}}_{k})V({\mathbf{x}}_{k})}{N}\right]$
$\displaystyle\displaystyle\leq(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})\overline{V}({\mathbf{x}}),$
(10)

where we used the fact that $\displaystyle\kappa^{\prime}$ and $\displaystyle
V$ are co-monotone and the FKG inequality for the last inequality.

For all $\displaystyle n\geq 1$, let $\displaystyle U_{n}=\\{{\mathbf{x}}\in
F,\ \max_{k\in\\{1,\ldots,|{\mathbf{x}}|\\}}{\mathbf{x}}_{k}>n\\}$ and define
the random variable $\displaystyle\theta_{n}:=\inf\\{t\geq 0,\
(X_{t}^{1},\ldots,X_{t}^{N})\in U_{n}\\}$. Note that $\displaystyle\theta_{n}$
is a stopping time since the constrained branching process is càdlàg.
Moreover, denoting by $\displaystyle C_{n}$ the set of points in
$\displaystyle F$ that can be reached in one jump from $\displaystyle
F\setminus U_{n}$, we observe that $\displaystyle C_{n}$ is bounded and hence
that $\displaystyle\bar{V}(X_{t\wedge\theta_{n}})$ is almost surely uniformly
bounded in $\displaystyle t\geq 0$. In particular, we deduce, using Dynkin’s
formula and the inequality (10), that, for all $\displaystyle n\geq 0$, all
$\displaystyle t\geq 0$ and all $\displaystyle{\mathbf{x}}\in F$,

$\displaystyle\mathbb{E}_{\mathbf{x}}(\bar{V}({\mathbf{X}}_{t\wedge\theta_{n}}))\leq\bar{V}({\mathbf{x}})+(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})\int_{0}^{t}\mathbb{E}_{x}(\bar{V}({\mathbf{X}}_{s\wedge\theta_{n}}))\,{\rm
d}s.$

Hence, using Grönwall’s inequality, we deduce that

$\displaystyle\displaystyle\mathbb{E}_{\mathbf{x}}(\bar{V}({\mathbf{X}}_{t\wedge\theta_{n}}))\leq\bar{V}({\mathbf{x}})e^{t(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})}.$

Since $\displaystyle\\{{\mathbf{X}}_{t\wedge\theta_{n}}\in
U_{n},\theta_{n}\leq t\\}=\\{\theta_{n}\leq t\\}$ up to a negligible event, we
conclude that

$\displaystyle\inf_{\mathbf{y}\in
U_{n}}\bar{V}(\mathbf{y})\mathbb{P}_{\mathbf{x}}(\theta_{n}\leq
t)\leq\bar{V}({\mathbf{x}})\,e^{t(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})}.$

Denote by $\displaystyle\tau_{\infty}$ the explosion time of the process.
Since the jump rate of the process is bounded up to time
$\displaystyle\theta_{n}$, we deduce that
$\displaystyle\theta_{n}\leq\tau_{\infty}$ almost surely. Hence

$\displaystyle\mathbb{P}_{\mathbf{x}}(\tau_{\infty}\leq
t)\leq\frac{\bar{V}({\mathbf{x}})e^{t(C+\|b\|_{\infty}+4\|\kappa^{\prime\prime}\|_{\infty})}}{\inf_{\mathbf{y}\in
U_{n}}\bar{V}(\mathbf{y})}\xrightarrow[n\to+\infty]{}0,$

so that $\displaystyle\tau_{\infty}=+\infty$ 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 $\displaystyle T>0$,
all bounded measurable functions $\displaystyle
f:\mathbb{Z}_{+}^{d}\to\mathbb{R}$ and all initial position
$\displaystyle{\mathbf{X}}(0)\in F$,

$\displaystyle\mathbb{E}_{{\mathbf{X}}(0)}\left[\Pi^{A}_{T}\Pi^{B}_{T}\sum_{k=1}^{|{\mathbf{X}}(t)|}f({\mathbf{X}}(t)_{k})\right]=\sum_{k=1}^{|{\mathbf{X}}(0)|}\mathbf{E}_{{\mathbf{X}}(0)_{k}}\left[f(X_{t})\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]$

where

$\displaystyle\Pi^{A}_{T}=\left(\frac{N_{min}-1}{N_{min}}\right)^{A_{T}}\text{
and }\Pi^{B}_{T}=\left(\frac{N_{max}+1}{N_{max}}\right)^{B_{T}}.$

In addition,

$\left\|\frac{\sum_{k=1}^{|{\mathbf{X}}(0)|}\mathbf{E}_{{\mathbf{X}}(0)_{k}}\left[f(X_{t})\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]}{\sum_{k=1}^{|{\mathbf{X}}(0)|}\mathbf{E}_{{\mathbf{X}}(0)_{k}}\left[\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]}-\frac{\sum_{k=1}^{|{\mathbf{X}}(t)|}f({\mathbf{X}}(t)_{k})}{|{\mathbf{X}}(t)|}\right\|_{2}\\\
\leq
C_{T}\frac{\|f\|_{\infty}}{\frac{1}{|{\mathbf{X}}(0)|}\sum_{k=1}^{|{\mathbf{X}}(0)|}\mathbf{E}_{{\mathbf{X}}(0)_{k}}\left[\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]}\,\frac{1}{\sqrt{|{\mathbf{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 $\displaystyle C_{1}>0,C_{2}>0,\eta_{0}>0,n_{0}\geq 0$
such that, for all $\displaystyle x\in\mathbb{Z}_{+}^{d}$ satisfying
$\displaystyle|x|\geq n_{0}$,

$\displaystyle\displaystyle\sum_{i=1}^{d}(\delta_{i}(x)-\beta_{i}(x))\geq
C_{1}\,|x|^{1+\eta_{0}}$ (11)
$\displaystyle\displaystyle\sum_{i=1}^{d}(\delta_{i}(x)-\beta_{i}(x))\geq
C_{2}(\|b\|_{\infty}-b(x)+\kappa(x))|x|.$ (12)

######  Proposition 17.

Assume that (11) and (12) hold true and that $\displaystyle X$ is irreducible.
Then there exists a probability distribution $\displaystyle\nu_{X}$ on
$\displaystyle\mathbb{Z}_{+}^{d}$ such that, for all probability measures
$\displaystyle\mu$ on $\displaystyle\mathbb{Z}_{+}^{d}$ and all bounded
measurable functions $\displaystyle f$ on $\displaystyle\mathbb{Z}_{+}^{d}$,

$\displaystyle\displaystyle\left|\frac{\mathbf{E}_{\mu}\left(f(X_{t})\exp\left(\int_{0}^{t}b(X_{s})-\kappa(X_{s})\,\mathrm{d}s\right)\right)}{\mathbf{E}_{\mu}\left(\exp\left(\int_{0}^{t}b(X_{s})-\kappa(X_{s})\,\mathrm{d}s\right)\right)}-\nu_{X}(f)\right|\leq
C_{X}\,e^{-\gamma_{X}t}\,\|f\|_{\infty},$ (13)

for some constants $\displaystyle C_{X}\geq 1$ and
$\displaystyle\gamma_{X}>0$. In addition, there exists a bounded function
$\displaystyle\eta_{X}:\mathbb{Z}_{+}^{d}\to(0,+\infty)$ and
$\displaystyle\lambda_{X}>-\|b\|_{\infty}$ such that for all probability
measures $\displaystyle\mu$ on $\displaystyle\mathbb{Z}_{+}^{d}$,

$\displaystyle\mathbf{E}_{\mu}\left(\eta_{X}(X_{t})\exp\left(\int_{0}^{t}b(X_{s})-\kappa(X_{s})\,\mathrm{d}s\right)\right)=e^{-\lambda_{X}t}\mu(\eta_{X}),\
\forall t\geq 0.$

###### Proof of Proposition 17.

First observe that, for all bounded measurable functions $\displaystyle
f:\mathbb{Z}_{+}^{d}\to\mathbb{R}$,

$\displaystyle\displaystyle\mathbf{E}_{\mu}\left(f(X_{t})\exp\left(\int_{0}^{t}b(X_{s})-\kappa(X_{s})\,\mathrm{d}s\right)\right)=e^{\|b\|_{\infty}t}\,\mathbf{E}^{b}_{\mu}\left(f(X^{b}_{t})\mathbf{1}_{t<\tau^{b}_{\partial}}\right),$

where $\displaystyle\mathbf{P}_{\mu}^{b}$ and
$\displaystyle\mathbf{E}_{\mu}^{b}$ denote the law and associated expectation
of the birth and death process $\displaystyle X^{b}$ with birth and death
parameters $\displaystyle\beta$ and $\displaystyle\delta$, killed at rate
$\displaystyle\kappa^{b}:=\|b\|_{\infty}-b+\kappa$ at a time denoted by
$\displaystyle\tau^{b}_{\partial}$. We now prove that there exist
$\displaystyle\nu_{X}$, $\displaystyle C_{X}>0$ and
$\displaystyle\gamma_{X}>0$ such that, for all bounded measurable functions
$\displaystyle f:\mathbb{Z}_{+}^{d}\to\mathbb{R}$

$\displaystyle\displaystyle\left|\mathbf{E}^{b}_{\mu}\left(\left.f(X^{b}_{t})\,\right|\,t<\tau^{b}_{\partial}\right)-\nu_{X}(f)\right|\leq
C_{x}\,e^{-\gamma_{X}t}\|f\|_{\infty},$ (14)

which is equivalent to (13).

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 $\displaystyle W(x)=|x|$ for all
$\displaystyle x\in\mathbb{Z}_{+}^{d}$. Assumption (11) implies that, for all
$\displaystyle x\in\mathbb{Z}_{+}^{p}$ large enough,

$\displaystyle\displaystyle\sum_{i=1}^{d}\left[\beta_{i}(x)(W(x+e_{i})-W(x))+\delta_{i}(x)(W(x-e_{i})-W(x))\right]\leq-
C_{1}W(x)^{1+\eta_{0}}.$

This classically entails that the birth and death process on
$\displaystyle\mathbb{Z}_{+}^{d}$ with birth rates $\displaystyle\beta$ and
death rates $\displaystyle\delta$ (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 [11]) and that Assumption 3 in [11] is
satisfied.

Moreover, choosing appropriately $\displaystyle a_{1}>1$, $\displaystyle
a_{2}>1$ and $\displaystyle\varepsilon>0$, and setting

$\displaystyle\displaystyle V(x)=\sum_{k=1}^{|x|}\frac{1}{k^{a_{1}}}\text{ and
}\varphi(x)=\sum_{k=|x|+1}^{+\infty}\frac{1}{k^{a_{2}}},$

one deduces, using the same computations as in Section 4 in [11], that
Assumption 1 of [11] is satisfied. The only difference is that the
infinitesimal generator that one needs to consider in the present case is

$\displaystyle\displaystyle
Lf(x)=\sum_{i=1}^{d}\left[\beta_{i}(x)(f(x+e_{i})-f(x))+\delta_{i}(x)(f(x-e_{i})-f(x))\right]-\kappa^{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 [11, Theorem 2.4]
applies to $\displaystyle X^{b}$, which concludes the proof of (13).

The last assertion of the proposition follows from [10, Proposition 2.3]: it
implies that there exists $\displaystyle\lambda_{0}>0$ and a bounded positive
function $\displaystyle\eta_{X}$ such that

$\displaystyle\mathbf{E}_{\mu}^{b}(\eta_{X}(X_{t}^{b}))=e^{-\lambda_{0}t}\eta_{X}.$

Multiplying by $\displaystyle e^{\|b\|_{\infty}t}$ and setting
$\displaystyle\lambda_{X}=\lambda_{0}-\|b\|_{\infty}$ concludes the proof. ∎

We consider now the situation where $\displaystyle\kappa$ is bounded and where
$\displaystyle N_{min}$ is large, and prove that the normalised empirical
distribution $\displaystyle\frac{m_{T}}{m_{T}(\mathbf{1}_{E})}$ of the
population health’s parameter converges toward $\displaystyle\nu_{X}$ when
$\displaystyle N_{min},N_{max}$ and $\displaystyle T$ go to infinity.

######  Proposition 18.

Assume that the conditions of Proposition 17 hold and that
$\displaystyle\kappa$ is bounded. Then, denoting by $\displaystyle m_{0}$ the
initial empirical measure of the population’s health parameter and setting
$\displaystyle\widetilde{m}_{0}=m_{0}/m_{0}(\mathbf{1}_{E})$, we have

$\sup_{T\in[0,+\infty)}\left\|\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}-\frac{\mathbf{E}_{\widetilde{m}_{0}}\left[f(X_{t})\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]}{\mathbf{E}_{\widetilde{m}_{0}}\left[\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right]}\right\|_{2}\leq(2C_{X}+C)N_{min}^{-\alpha/2}\,\|f\|_{\infty},$

where
$\displaystyle{\displaystyle\alpha=\gamma_{X}/((c+1)\|b\|_{\infty}+\|\kappa\|_{\infty}+\gamma_{X})\in(0,1)}$,
$\displaystyle C_{X},\gamma_{X}$ are from Proposition 17, and $\displaystyle
C,c$ are from Theorem 6. In particular,

$\displaystyle\left\|\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}-\nu_{X}(f)\right\|_{2}\leq\left[(2C_{X}+C)N_{min}^{-\alpha/2}\,\|f\|_{\infty}+C_{X}\,e^{-\gamma_{X}T}\right]\,\|f\|_{\infty}.$

######  Remark 19.

This result can be adapted to include the case where $\displaystyle 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 [23, 10, 9, 36, 25]),
provided one can prove that the particle system does not degenerate toward the
boundary. However this can be challenging, see for instance [33, 56, 5].

###### Proof of Proposition 18.

We assume without loss of generality that $\displaystyle\|f\|_{\infty}=1$ and
that $\displaystyle N_{min}>\sqrt{N_{max}}$ (otherwise the result is trivial).
We define the semigroup $\displaystyle Q$ by

$\displaystyle
Q_{t}f(x)=\mathbf{E}_{x}\left[f(X_{t})\exp\left(\int_{0}^{t}(b(X_{s})-\kappa(X_{s}))\,\mathrm{d}s\right)\right].$

We have, by Proposition 17 applied first to
$\displaystyle\mu=m_{0}Q_{T-t}/m_{0}Q_{T-t}\mathbf{1}_{E}$ and then to
$\displaystyle\widetilde{m}_{T-t}=m_{T-t}/m_{T-t}\mathbf{1}_{E}$, and by
Theorem 6, for all $\displaystyle 0\leq t\leq T$,

$\displaystyle\displaystyle\left\|\frac{m_{0}Q_{T}f}{m_{0}Q_{T}\mathbf{1}_{E}}-\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}\right\|_{2}$
$\displaystyle\displaystyle\leq\left\|\frac{m_{0}Q_{T-t}Q_{t}f}{m_{0}Q_{T-t}Q_{t}\mathbf{1}_{E}}-\frac{m_{T-t}(Q_{t}f)}{m_{T-t}(Q_{t}\mathbf{1}_{E})}\right\|_{2}+\left\|\frac{m_{T-t}(Q_{t}f)}{m_{T-t}(Q_{t}\mathbf{1}_{E})}-\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}\right\|_{2}$
$\displaystyle\displaystyle\leq
2C_{X}e^{-\gamma_{X}t}+C\exp(c\|b\|_{\infty}t)\left\|\frac{\|Q_{t}\mathbf{1}_{E}\|_{\infty}}{\widetilde{m}_{T-t}Q_{t}\mathbf{1}_{E}}\frac{1}{\sqrt{m_{T-t}\mathbf{1}_{E}}}\right\|_{2},$
(15)

where

$\displaystyle\displaystyle\left\|\frac{\|Q_{t}\mathbf{1}_{E}\|_{\infty}}{\widetilde{m}_{T-t}Q_{t}\mathbf{1}_{E}}\frac{1}{\sqrt{m_{T-t}\mathbf{1}_{E}}}\right\|_{2}$
$\displaystyle\displaystyle\leq\frac{\exp(\|b\|_{\infty}t)}{\exp(-\|\kappa\|_{\infty}t)}\frac{1}{\sqrt{N_{min}}}.$

Set

$\displaystyle
t_{0}=\frac{\ln\left(\sqrt{N_{min}}\right)}{(c+1)\|b\|_{\infty}+\|\kappa\|_{\infty}+\gamma_{X}},$

so that $\displaystyle
e^{-\gamma_{X}t_{0}}=e^{((c+1)\|b\|_{\infty}+\|\kappa\|_{\infty})t_{0}}/\sqrt{N_{min}}$.
If $\displaystyle T\geq t_{0}$, then choosing $\displaystyle t=t_{0}$ in (15),
we get

$\displaystyle\displaystyle\left\|\frac{m_{0}Q_{T}f}{m_{0}Q_{T}\mathbf{1}_{E}}-\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}\right\|_{2}$
$\displaystyle\displaystyle\leq(2C_{X}+C)e^{-\gamma_{X}t_{0}}$
$\displaystyle\displaystyle=(2C_{X}+C)N_{min}^{-\alpha/2},$

with
$\displaystyle\alpha=\frac{\gamma_{X}}{(c+1)\|b\|_{\infty}+\|\kappa\|_{\infty}+\gamma_{X}}$.
For $\displaystyle T\leq t_{0}$, Theorem 6 implies that

$\displaystyle\displaystyle\left\|\frac{m_{0}Q_{T}f}{m_{0}Q_{T}\mathbf{1}_{E}}-\frac{m_{T}(f)}{m_{T}(\mathbf{1}_{E})}\right\|_{2}$
$\displaystyle\displaystyle\leq
C\exp\left((c+1)\|b\|_{\infty}t_{0}+\|\kappa\|_{\infty}t_{0}\right)\frac{1}{\sqrt{N_{min}}}$
$\displaystyle\displaystyle=CN_{min}^{-\alpha/2}$

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. ∎

###  3.2 Piecewise deterministic Markov processes

Due to Assumption 1 in Section 2, the results presented thus far do not hold
for the case where $\displaystyle(X_{t})_{t\geq 0}$ is a piecewise
deterministic Markov process (PDMP) and $\displaystyle E$ is a bounded domain
with absorbing boundary conditions, since two independent copies of the
process may hit $\displaystyle\partial 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 $\displaystyle E=D\times V$ where $\displaystyle
D\subset\mathbb{R}^{3}$ is the spatial domain, which we assume to be open and
bounded with smooth boundary, $\displaystyle\partial D$, and $\displaystyle
V=\\{\upsilon\in\mathbb{R}^{3}:v_{min}<|\upsilon|<v_{max}\\}$ is the velocity
domain, where $\displaystyle 0<v_{min}\leq v_{max}<\infty$. The NTE can then
be stated as follows:

$\displaystyle\displaystyle\frac{\partial}{\partial t}\psi_{t}(r,\upsilon)$
$\displaystyle\displaystyle=\upsilon\cdot\nabla\psi_{t}(r,\upsilon)-(\sigma_{\mathtt{s}}(r,\upsilon)+\sigma_{\mathtt{f}}(r,\upsilon))\psi_{t}(r,\upsilon)$
$\displaystyle\displaystyle+\sigma_{\texttt{s}}(r,\upsilon)\int_{V}\psi_{t}(r,\upsilon^{\prime})\pi_{\texttt{s}}(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime}+\sigma_{\texttt{f}}(r,\upsilon)\int_{V}\psi_{t}(r,\upsilon^{\prime})\pi_{\texttt{f}}(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime},$ (16)

where

$\displaystyle\displaystyle\sigma_{\texttt{s}}(r,\upsilon)$
$\displaystyle\displaystyle:\text{ the rate at which scattering occurs from
incoming velocity $\displaystyle\upsilon$ at position $\displaystyle r$,}$
$\displaystyle\displaystyle\sigma_{\texttt{f}}(r,\upsilon)$
$\displaystyle\displaystyle:\text{ the rate at which fission occurs from
incoming velocity $\displaystyle\upsilon$ at position $\displaystyle r$,}$
$\displaystyle\displaystyle\pi_{\texttt{s}}(r,\upsilon,\upsilon^{\prime})$
$\displaystyle\displaystyle:\text{ probability density that an incoming
velocity $\displaystyle\upsilon$ at position $\displaystyle r$ scatters to an
}$ $\displaystyle\displaystyle\hskip 14.22636pt\text{outgoing velocity, with
probability $\displaystyle\upsilon^{\prime}$ satisfying
}\textstyle{\int_{V}}\pi_{\texttt{s}}(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime}=1,\text{ and }$
$\displaystyle\displaystyle\pi_{\texttt{f}}(r,\upsilon,\upsilon^{\prime})$
$\displaystyle\displaystyle:\text{ density of expected neutron yield at
velocity $\displaystyle\upsilon^{\prime}$ from fission with }$
$\displaystyle\displaystyle\hskip 14.22636pt\text{incoming velocity
$\displaystyle\upsilon$ satisfying
}\textstyle{\int_{V}\pi_{\texttt{f}}}(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime}<\infty.$

Moreover, we impose the following initial and boundary conditions

$\left.\begin{array}[]{ll}\psi_{0}(r,\upsilon)=g(r,\upsilon)&\text{ for }r\in
D,\upsilon\in{V},\\\ &\\\ \psi_{t}(r,\upsilon)=0&\text{ for }r\in\partial
D\text{ if }\upsilon\cdot{\bf n}_{r}>0,\end{array}\right\\}$ (17)

where $\displaystyle{\bf n}_{r}$ is the outward facing normal of
$\displaystyle D$ at $\displaystyle r\in\partial D$ and $\displaystyle
g:D\times{V}\to[0,\infty)$ is a bounded, measurable function.

Under the assumptions,

  * (A1)

$\displaystyle\sigma_{\mathtt{s}}$, $\displaystyle\sigma_{\mathtt{f}}$,
$\displaystyle\pi_{\mathtt{s}}$ and $\displaystyle\pi_{\mathtt{f}}$ are
uniformly bounded away from infinity.

  * (A2)

$\displaystyle\inf_{r\in D,\upsilon,\upsilon^{\prime}\in
V}(\sigma_{\texttt{s}}(r,\upsilon)\pi_{\texttt{s}}(r,\upsilon,\upsilon^{\prime})+\sigma_{\texttt{f}}(r,\upsilon)\pi_{\texttt{f}}(r,\upsilon,\upsilon^{\prime}))>0$,

it was shown in [14, 36] that the NTE can be modelled using a weighted PDMP,
also known as the neutron random walk (NRW). More precisely, setting

$\displaystyle\displaystyle\alpha(r,\upsilon)$
$\displaystyle\displaystyle=\sigma_{\texttt{s}}(r,\upsilon)+\sigma_{\texttt{f}}(r,\upsilon)\int_{V}\pi_{\texttt{f}}(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime}$ (18)
$\displaystyle\displaystyle\pi(r,\upsilon,\upsilon^{\prime})$
$\displaystyle\displaystyle=(\alpha(r,\upsilon))^{-1}\left[\sigma_{\texttt{s}}(r,\upsilon)\pi_{\texttt{s}}(r,\upsilon,\upsilon^{\prime})+\sigma_{\texttt{f}}(r,\upsilon)\pi_{\texttt{f}}(r,\upsilon,\upsilon^{\prime})\right],$
(19) $\displaystyle\displaystyle\beta(r,\upsilon)$
$\displaystyle\displaystyle=\alpha(r,\upsilon)-\sigma_{\texttt{s}}(r,\upsilon)-\sigma_{\texttt{f}}(r,\upsilon),$
(20)

we define the NRW $\displaystyle((R_{t},\Upsilon_{t})_{t\geq
0},\mathbf{P}_{(r,\upsilon)})$ as follows. From an initial configuration
$\displaystyle(r,\upsilon)\in D\times V$, the particle will propagate linearly
until it is either absorbed at the boundary of the domain or, at rate
$\displaystyle\alpha$ the process scatters off a nucleus and a new velocity is
chosen according to $\displaystyle\pi$. Then, for a bounded, measurable
function $\displaystyle g$, and $\displaystyle(r,\upsilon)\in D\times V$

$\displaystyle\psi_{t}[g](r,\upsilon)\coloneqq{\mathbf{E}}_{(r,\upsilon)}\left[{\rm
e}^{\int_{0}^{t}\beta(R_{s},\Upsilon_{s}){\rm
d}s}g(R_{t},\Upsilon_{t})\mathbf{1}_{(t<\tau_{D})}\right],$

solves (16)–(17), where $\displaystyle\tau_{D}$ is the first exit time of the
NRW from $\displaystyle D$.

It was also shown in [36] that if (A1) and (A2) are satisfied, and that if
there exists an $\displaystyle\varepsilon>0$ such that

  * (B1)

$\displaystyle\textstyle{D_{\varepsilon}\coloneqq\\{r\in D:\inf_{y\in\partial
D}|r-y|>\varepsilon{\texttt{v}}_{\texttt{max}}\\}}$ is non-empty and
connected,

  * (B2)

there exist $\displaystyle 0<s_{\varepsilon}<t_{\varepsilon}$ and
$\displaystyle{\color[rgb]{0,0,0}\gamma}>0$ such that, for all $\displaystyle
r\in D\backslash D_{\varepsilon}$, there exists $\displaystyle K_{r}\subset V$
measurable such that
$\displaystyle\text{Vol}({K_{r}})\geq{\color[rgb]{0,0,0}\gamma}>0$ and for all
$\displaystyle\upsilon\in K_{r}$, $\displaystyle r+\upsilon s\in
D_{\varepsilon}$ for every $\displaystyle
s\in[s_{\varepsilon},t_{\varepsilon}]$ and $\displaystyle r+\upsilon
s\notin\partial D$ for all $\displaystyle s\in[0,s_{\varepsilon}]$,

the semigroup $\displaystyle(\psi_{t})_{t\geq 0}$ exhibits the following
Perron Frobenius behaviour.

######  Theorem ([36]).

There exists a $\displaystyle\lambda_{*}\in\mathbb{R}$, a positive right
eigenfunction $\displaystyle\varphi\in L^{+}_{\infty}(D\times V)$ and a left
eigenmeasure which is absolutely continuous with respect to Lebesgue measure
on $\displaystyle D\times V$ with density $\displaystyle\tilde{\varphi}\in
L^{+}_{\infty}(D\times V)$, both having associated eigenvalue
$\displaystyle{\rm e}^{\lambda_{*}t}$, and such that $\displaystyle\varphi$
(resp. $\displaystyle\tilde{\varphi}$) is uniformly (resp. a.e. uniformly)
bounded away from zero on each compactly embedded subset of $\displaystyle
D\times V$. In particular, for all $\displaystyle g\in L^{+}_{\infty}(D\times
V)$,

$\langle\tilde{\varphi},\psi_{t}[g]\rangle={\rm
e}^{\lambda_{*}t}\langle\tilde{\varphi},g\rangle\quad\text{(resp.
}\psi_{t}[\varphi]={\rm e}^{\lambda_{*}t}\varphi\text{)},\quad t\geq 0.$ (21)

Moreover, there exists $\displaystyle\varepsilon>0$ such that

$\sup_{g\in L^{+}_{\infty}(D\times V):\|g\|_{\infty}\leq 1}\left\|{\rm
e}^{-\lambda_{*}t}{\varphi}^{-1}{\psi_{t}[g]}-\langle\tilde{\varphi},g\rangle\right\|_{\infty}=O({\rm
e}^{-\varepsilon t}),\qquad t\geq 0.$ (22)

This theorem therefore provides a way to build Monte Carlo simulations to
estimate the eigenvalue $\displaystyle\lambda_{*}$ and the eigenfunctions
$\displaystyle\varphi$ and $\displaystyle\tilde{\varphi}$ via the NRW. We
refer the reader to [15] 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 ‘$\displaystyle h$-transform’ was developed in [15],
where one biases the NRW to prevent it from exiting the domain.

More precisely, let $\displaystyle h$ be a bounded positive function on
$\displaystyle D\times V$ such that $\displaystyle h=0$ on
$\displaystyle\partial D$ and define the following change of measure,

$\left.\frac{{\rm d}\mathbf{P}^{h}_{(r,\upsilon)}}{{\rm
d}\mathbf{P}_{(r,\upsilon)}}\right|_{\sigma((R_{s},\Upsilon_{s}),s\leq
t)}:=\exp\left(-\int_{0}^{t}\frac{{\emph{{J}}}h(R_{s},\Upsilon_{s})}{h(R_{s},\Upsilon_{s})}{\rm
d}s\right)\prod_{i=1}^{N_{t}}\frac{h(R_{T_{i}},\Upsilon_{T_{i}})}{h(R_{T_{i}},\Upsilon_{T_{i-1}})}\,\mathbf{1}_{(t<\tau^{D})},$
(23)

where

${\emph{{J}}}g(r,\upsilon)=\alpha(r,\upsilon)\int_{V}[g(r,\upsilon^{\prime})-g(r,\upsilon)]\pi(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime},$ (24)

$\displaystyle(T_{i})_{i\geq 0}$ are the scatter times of the NRW with
$\displaystyle T_{0}=0$, and $\displaystyle N_{t}=\sup\\{i:T_{i}\leq t\\}$.

Then, $\displaystyle((R,\Upsilon),\mathbf{P}^{h})$ also defines a NRW but with
scattering operator

${\emph{{J}}}_{h}g(r,\upsilon)=\alpha(r,\upsilon)\int_{V}[g(r,\upsilon^{\prime})-g(r,\upsilon)]\frac{h(r,\upsilon^{\prime})}{h(r,\upsilon)}\pi(r,\upsilon,\upsilon^{\prime}){\rm
d}\upsilon^{\prime}.$ (25)

The idea is that the factor $\displaystyle
h(r,\upsilon^{\prime})/h(r,\upsilon)$ 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 $\displaystyle h$-transformed
process to estimate the leading eigentriple
$\displaystyle(\lambda_{*},\varphi,\tilde{\varphi})$, for example. Indeed, it
is straightforward to show that the leading eigentriple of the $\displaystyle
h$-transformed process is given by
$\displaystyle(\lambda_{*},\varphi/h,h\tilde{\varphi})$ and so we may use the
BBMMI process to estimate this triple, and hence obtain an estimate for the
original quantities $\displaystyle(\lambda_{*},\varphi,\tilde{\varphi})$. 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 $\displaystyle\delta=1/(2\|\alpha\|_{\infty})$ and let
$\displaystyle\phi:[0,+\infty)\to[0,\delta]$ be a regular non-decreasing
function such that $\displaystyle\phi(x)=x$ for all $\displaystyle
x\in[0,\delta/2]$ and $\displaystyle\phi(x)=\delta$ for all $\displaystyle
x\geq\delta$. We set

$\displaystyle h(r,\upsilon)=\phi\left(\kappa_{r,\upsilon}\right).$

Using the following result, one can express expectations with respect to
$\displaystyle\mathbf{P}$ as Feynman-Kac formulas related to
$\displaystyle\mathbf{P}^{h}$. In what follows, we say that
$\displaystyle(r,\upsilon)$ is in a vicinity of $\displaystyle\partial D$ if
$\displaystyle h(r,\upsilon)<\delta/2$.

######  Lemma 20.

Suppose (A1), (A2) hold. A process with law
$\displaystyle\mathbf{P}^{h}_{(r,v)}$ does not hit the boundary with
probability one. Moreover, setting $\displaystyle
b=\left(\frac{Lh}{h}\right)_{+}$ and
$\displaystyle\kappa=\left(\frac{Lh}{h}\right)_{-}$, the function
$\displaystyle b$ is bounded on $\displaystyle D\times V$ and
$\displaystyle\kappa$ is locally bounded. Let
$\displaystyle\mathbf{P}^{\kappa}_{r,\upsilon}$ be the law of a process with
law $\displaystyle\mathbf{P}^{h}_{r,\upsilon}$ with additional soft killing at
rate $\displaystyle\kappa$, then we have

$\displaystyle\displaystyle\mathbf{E}_{(r,\upsilon)}\left[f(R_{t},\Upsilon_{t})\,\mathbf{1}_{(t<\tau^{D})}\right]=h(r,\upsilon)\,\mathbf{E}^{\kappa}_{(r,\upsilon)}\left[\exp\left(\int_{0}^{t}b(R_{s},\Upsilon_{s})\,{\rm
d}s\right)\frac{f(R_{t},\Upsilon_{t})}{h(R_{t},\Upsilon_{t})}\mathbf{1}_{(t<\tau_{\kappa})}\right],$
(26)

where $\displaystyle\tau_{\kappa}$ denotes the soft killing time of the
process with law $\displaystyle\mathbf{P}^{\kappa}_{r,\upsilon}$.

###### Proof.

We have

$\displaystyle\displaystyle\psi_{t}[g](r,\upsilon)$
$\displaystyle\displaystyle=\mathbf{E}^{h}_{(r,\upsilon)}\left[\exp\left(\int_{0}^{t}\frac{{\emph{{J}}}h(R_{s},\Upsilon_{s})}{h(R_{s},\Upsilon_{s})}+\beta(R_{s},\Upsilon_{s}){\rm
d}s\right)\prod_{i=1}^{N_{t}}\frac{h(R_{T_{i}},\Upsilon_{T_{i-1}})}{h(R_{T_{i}},\Upsilon_{T_{i}})}g(R_{t},\Upsilon_{t})\mathbf{1}_{(t<\tau^{D})}\right]$
(27)

Note that

$\inf_{r\in D,\upsilon\in
V}\lim_{s\to\kappa_{r,\upsilon}^{D}}\frac{|\upsilon|(\kappa_{r,\upsilon}^{D}-s)}{h(r+\upsilon
s,\upsilon)}>0,$ (28)

where $\displaystyle\kappa_{r,\upsilon}^{D}:=\inf\\{t>0:r+\upsilon t\not\in
D\\}$. Hence, Theorem 7.1 in [15] implies that, under
$\displaystyle\mathbf{P}^{h}_{(r,\upsilon)}$, $\displaystyle r\in
D,\upsilon\in V$, $\displaystyle(R,\Upsilon)$ does not hit the boundary with
probability one.

Recall that the infinitesimal generator associated to
$\displaystyle((R_{t},\Upsilon_{t})_{t\geq 0},\mathbf{P})$ is given by

$\displaystyle
Lf(r,\upsilon)=\mathtt{T}f+\alpha(r,\upsilon)\int_{V}(f(r,\upsilon^{\prime})-f(r,\upsilon))\,\pi(r,\upsilon,\upsilon^{\prime}),{\rm
d}\upsilon,$

where $\displaystyle\mathtt{T}$ is defined for regular functions by
$\displaystyle\mathtt{T}f(r,\upsilon)=\upsilon\cdot\nabla_{r}f(r,\upsilon)$
and, more generally and when the limit exists, by

$\displaystyle\mathtt{T}f(r,\upsilon)=\lim_{\varepsilon\to
0^{+}}\frac{f(r+\varepsilon\upsilon)-f(r,\upsilon)}{\varepsilon}.$

In particular, the conditions of Remark 7.1 of [15] are satisfied and hence

$\left.\frac{{\rm d}\mathbf{P}^{h}_{(r,\upsilon)}}{{\rm
d}\mathbf{P}_{(r,\upsilon)}}\right|_{\sigma((R_{s},\Upsilon_{s}),s\leq
t)}=\exp\left(-\int_{0}^{t}\frac{Lh(R_{s},\Upsilon_{s})}{h(R_{s},\Upsilon_{s})}{\rm
d}s\right)\frac{h(R_{t},\Upsilon_{t})}{h(r,\upsilon)}\,\mathbf{1}_{(t<\tau^{D})}.$
(29)

Finally, we obtain, for all measurable functions $\displaystyle f:D\times
V\to[0,+\infty)$,

$\displaystyle\displaystyle\mathbf{E}_{(r,\upsilon)}\left[f(R_{t},\Upsilon_{t})\,\mathbf{1}_{(t<\tau^{D})}\right]=h(r,\upsilon)\,\mathbf{E}^{h}_{(r,\upsilon)}\left[\exp\left(\int_{0}^{t}\frac{Lh(R_{s},\Upsilon_{s})}{h(R_{s},\Upsilon_{s})}{\rm
d}s\right)\frac{f(R_{t},\Upsilon_{t})}{h(R_{t},\Upsilon_{t})}\right]$ (30)

where have used the fact that $\displaystyle\mathbf{1}_{(t<\tau^{D})}=1$,
$\displaystyle\mathbf{P}^{h}_{(r,\upsilon)}$-almost surely, deduced from (28).
Now, introducing $\displaystyle b=\left(\frac{Lh}{h}\right)_{+}$ and
$\displaystyle\kappa=\left(\frac{Lh}{h}\right)_{-}$, and considering the law
$\displaystyle\mathbf{P}^{\kappa}_{r,\upsilon}$ of a process with law
$\displaystyle\mathbf{P}^{h}_{r,\upsilon}$ with soft killing at rate
$\displaystyle\kappa$, we deduce that

$\displaystyle\displaystyle\mathbf{E}_{(r,\upsilon)}\left[f(R_{t},\Upsilon_{t})\,\mathbf{1}_{(t<\tau^{D})}\right]=h(r,\upsilon)\,\mathbf{E}^{\kappa}_{(r,\upsilon)}\left[\exp\left(\int_{0}^{t}b(R_{s},\Upsilon_{s})\,{\rm
d}s\right)\frac{f(R_{t},\Upsilon_{t})}{h(R_{t},\Upsilon_{t})}\mathbf{1}_{(t<\tau_{\kappa})}\right],$
(31)

where $\displaystyle\tau_{\kappa}$ denotes the soft killing time of the
process (with law $\displaystyle\mathbf{P}^{h}_{r,\upsilon}$ with soft killing
at rate $\displaystyle\kappa$).

Now, we observe that, in a vicinity of $\displaystyle\partial D$,

$\displaystyle\displaystyle Lh(r,\upsilon)$
$\displaystyle\displaystyle=-1+\alpha(r,\upsilon)\int_{V}(h(r,\upsilon^{\prime})-h(r,\upsilon))\,\pi(r,\upsilon,\upsilon^{\prime})\,{\rm
d}\upsilon^{\prime}$ (32)
$\displaystyle\displaystyle\leq-1+\|\alpha\|_{\infty}\delta=-1/2$

and hence that

$\displaystyle\sup_{(r,\upsilon)\in D\times
V}\frac{Lh(r,\upsilon)}{h(r,\upsilon)}<+\infty.$

This implies that $\displaystyle b$ is bounded. The fact that
$\displaystyle\kappa$ is locally bounded is an immediate consequence of the
regularity of $\displaystyle 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
$\displaystyle\mathbf{P}$. Assumption 1 is clearly satisfied by a process with
law $\displaystyle\mathbf{P}^{\kappa}_{(r,\upsilon)}$, 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
$\displaystyle p^{i}$ results in a size constrained process, i.e. there is a
constant $\displaystyle N_{max}$ such that $\displaystyle\bar{N}_{t}\leq
N_{max}$ for all $\displaystyle t\geq 0$.

######  Proposition 21.

Under the assumptions (A1), (A2), (B1) and (B2), any  bounded BBMMI driven by
$\displaystyle\mathbf{P}^{\kappa}_{(r,\upsilon)}$ satisfies Assumption 2.

###### Proof.

We denote by $\displaystyle\sigma_{1}<\sigma_{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
$\displaystyle\sigma_{n}$ is well defined since, almost surely, the number of
event occurring before a time $\displaystyle T$ goes to infinity when
$\displaystyle T\to+\infty$. 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 $\displaystyle 0$.

From the expression of the scattering operator $\displaystyle{\emph{{J}}}_{h}$
in (25), one observes that the scattering rate toward a direction in a set
$\displaystyle V^{\prime}\subset V$ is given by
$\displaystyle\int_{V^{\prime}}\alpha(r,\upsilon)\frac{h(r,\upsilon^{\prime})}{h(r,\upsilon)}\pi(r,\upsilon,\upsilon^{\prime})\mathrm{d}\upsilon^{\prime}$.
But, under the assumptions (B1–2), there exists $\displaystyle\varepsilon>0$,
$\displaystyle\beta>0$ and $\displaystyle\eta>0$ such that, for any point
$\displaystyle r\in D$, there exists $\displaystyle V_{r}$ with $\displaystyle
r+[0,\beta]V_{r}\subset D$ and
$\displaystyle\text{Leb}_{d}((r+[0,\beta]V_{r})\cap\\{r^{\prime}\in D,\
d(r^{\prime},\partial D)\geq\varepsilon\\})>\eta$, where
$\displaystyle\text{Leb}_{d}$ denotes the $\displaystyle d$-dimensional
Lebesgue measure. In particular, using the fact that $\displaystyle
h(r,\upsilon^{\prime})$ is lower bounded on $\displaystyle V_{r}$ (uniformly
over $\displaystyle r\in D$). We deduce that there exists a constant
$\displaystyle\underline{s}>0$ such that the scattering rate toward a
direction in $\displaystyle V_{r}$ is lower bounded by
$\displaystyle\underline{s}/h(r,\upsilon)$, uniformly in
$\displaystyle(r,\upsilon)\in D\times V$.

Using (32) to bound $\displaystyle Lh$ from below, there exists a constant
$\displaystyle\bar{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 $\displaystyle\bar{c}/h(r,\upsilon)$, uniformly in
$\displaystyle(r,\upsilon)\in D\times V$. Hence each event has a probability
greater than $\displaystyle\underline{s}/\bar{c}$ to be a scattering event
toward a direction in $\displaystyle V_{r}$, independently of the past of the
process, and in particular independently of the time event
$\displaystyle\sigma_{n}$: formally, for all $\displaystyle n\geq 0$,

$\displaystyle\mathbf{P}^{\kappa}\left(\sigma_{n}\text{ is a scattering event
and }(R_{\sigma_{n}},\Upsilon_{\sigma_{n}})\in\\{R_{\sigma_{n}-}\\}\times
V_{R_{\sigma_{n}-}}\mid\sigma_{n},R_{\sigma_{n}-},\Upsilon_{\sigma_{n}-}\right)\geq\underline{s}\,/\bar{c}>0,$

where $\displaystyle(R,\Upsilon)$ 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 $\displaystyle\\{r^{\prime}\in D,\ d(r^{\prime},\partial
D)\geq\varepsilon\\}$ and the fact that the total branching rate of the system
is uniformly bounded, we deduce that there exists a constant
$\displaystyle\underline{p}>0$ and a time $\displaystyle T>0$ such that, for
all $\displaystyle r\in D$ and $\displaystyle\upsilon\in V_{r}$,
$\displaystyle{\mathbb{P}}(\sigma_{1}\geq T)\geq\underline{p}$.

Hence, we deduce from the strong Markov property that

$\displaystyle\displaystyle\mathbf{P}^{\kappa}(\sigma_{n+1}<T)$
$\displaystyle\displaystyle\leq\mathbf{E}^{\kappa}\left[{\mathbf{1}}_{\sigma_{n}<T}\,\mathbf{P}^{\kappa}_{X_{\sigma_{n}}}(\sigma_{1}<T)\right]$
$\displaystyle\displaystyle\leq\mathbf{E}^{\kappa}\left[{\mathbf{1}}_{\sigma_{n}<T}\,\left(1-{\mathbf{1}}_{\Upsilon_{\sigma_{n}}\in
V_{R_{\sigma_{n}}}}\mathbf{P}^{\kappa}_{X_{\sigma_{n}}}(\sigma_{1}<T)\right)\right]$
$\displaystyle\displaystyle\leq\mathbf{E}^{\kappa}\left[{\mathbf{1}}_{\sigma_{n}<T}\,\left(1-{\mathbf{1}}_{\Upsilon_{\sigma_{n}}\in
V_{R_{\sigma_{n}}}}\underline{p}\right)\right]$
$\displaystyle\displaystyle\leq\mathbf{P}^{\kappa}(\sigma_{n}<T)(1-\underline{s}\,/\bar{c}\,\underline{p}).$

This shows that the number of events occurring before time $\displaystyle T$
is stochastically dominated by a geometric random variable with parameter
$\displaystyle\underline{s}\,/\bar{c}\,\underline{p}>0$, so that it is finite
almost surely and that its expectation is bounded by a constant
$\displaystyle\bar{a}$. Since this is true for all particles, we deduce that
total number of events has an expectation bounded by $\displaystyle
N_{max}\bar{a}$. Finally, observing that this bound does not depend on the
initial position, we deduce using the Markov property at times $\displaystyle
T,2T,\ldots$ that the number of event occurring before a time horizon
$\displaystyle nT$ is bounded by $\displaystyle kN_{max}\bar{a}$, and hence
that it is finite almost surely.

This concludes the proof of Proposition 21. ∎

###  3.3 Numerical properties of the constrained branching process

The aim of this section is to discuss some of the numerical properties of the
$\displaystyle N_{min}$–$\displaystyle 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 [18, 19]. 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
$\displaystyle\Pi_{T}^{A}\Pi^{B}_{T}$ has high variance.

####  3.3.1 Comparison with the FV IPS

We recall that the FV IPS simulates a set of $\displaystyle 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
$\displaystyle N-1$ particles and duplicating it. Note that the FV IPS is a
particular instance of the $\displaystyle N_{min}$–$\displaystyle N_{max}$
IPS, with a null branching rate and $\displaystyle N_{min}=N_{max}=N$. In
particular, Theorem 6 applies to the FV IPS (see [57] 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 $\displaystyle X=(X_{t})_{t\geq 0}$, under fairly weak
assumptions, the following many-to-one lemma holds

$\displaystyle\psi_{t}[f](x)\coloneqq\mathbb{E}_{\delta_{x}}[\langle
f,X_{t}\rangle]=\mathbf{E}_{x}\left[{\rm e}^{\int_{0}^{t}\beta(Y_{s}){\rm
d}s}f(Y_{t})\right],$

where $\displaystyle\beta$ is determined by the branching rate and mean
offspring of $\displaystyle X$, and the law of $\displaystyle Y$ is uniform
amongst all paths of $\displaystyle X$. In the case where $\displaystyle\beta$
is uniformly bounded above, one may compensate the right-hand side above by
its maximum, $\displaystyle\bar{\beta}$, to introduce a penalised process:
$\displaystyle\psi_{t}^{\mathtt{k}}[f](x)\coloneqq{\rm
e}^{-\bar{\beta}t}\psi_{t}[f](x)$. The penalisation has the effect of killing
$\displaystyle Y$ at an additional rate $\displaystyle\bar{\beta}-\beta$, thus
yielding a sub-Markov process with semigroup
$\displaystyle\psi_{t}^{\mathtt{k}}$. This then allows one to use the FV IPS
to obtain estimates of this semigroup and, in some situations, of its
associated quasi-stationary 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 [6, 18, 51, 32, 56,
34, 8].

In this section, we show on a simple example that the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ algorithm allows one to overcome some of the
limitations of the FV IPS. For simplicity, we only consider the situation
where $\displaystyle N_{min}=N_{max}$.

##### A simple branching birth and death process.

We consider a birth and death process on $\displaystyle E=\\{1,\dots,M\\}$ for
some parameter $\displaystyle M\in\mathbb{N}$. For a particle occupying state
$\displaystyle x$, one of the following things may occur:

  * •

the particle jumps with rate $\displaystyle x^{2}$, to state $\displaystyle
1\vee(x-1)$ with probability $\displaystyle\frac{x}{x+1}$ and to state
$\displaystyle(x+1)\wedge M$ with probability $\displaystyle\frac{1}{x+1}$,

  * •

at rate $\displaystyle 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 $\displaystyle(\psi_{t})_{t\geq 0}$ admits a limiting
distribution up to an exponential rescaling (see for instance [16]): there
exist a constant $\displaystyle\lambda_{M}\in\mathbb{R}$, a positive function
$\displaystyle\eta_{M}:\\{1,\ldots,M\\}\to(0,+\infty)$ and a probability
measure $\displaystyle\nu_{M}$ on $\displaystyle\\{1,\ldots,M\\}$ such that

$\displaystyle\displaystyle\sup_{f:\\{1,\ldots,M\\},\|f\|_{\infty}\leq
1}\left|e^{-\lambda_{M}t}\psi_{t}[f](x)-\eta_{M}(x)\nu_{M}(f)\right|\leq
C_{M}e^{-\gamma_{M}t},\qquad t\geq 0,$

for some positive constants $\displaystyle C_{M},\gamma_{M}$.

The associated birth and death process with killing evolves on
$\displaystyle\\{1,\ldots,M\\}$ with the same jump rates, without branching,
and with the additional killing rate $\displaystyle\kappa(x)=M-x$. Denoting as
above by $\displaystyle\psi^{\mathtt{k}}$ the associated sub-Markov semigroup,
we have $\displaystyle\psi^{\mathtt{k}}_{t}=e^{-Mt}\psi_{t}$ and hence

$\displaystyle\displaystyle\sup_{f:\\{1,\ldots,M\\},\|f\|_{\infty}\leq
1}\left|e^{-(\lambda_{M}-M)t}\psi^{\mathtt{k}}_{t}[f](x)-\eta_{M}(x)\nu_{M}(f)\right|\leq
C_{M}e^{-\gamma_{M}t}.$

##### Theoretical convergence of the $\displaystyle N_{min}$–$\displaystyle
N_{max}$ and of the FV IPSs.

In what follows, $\displaystyle m^{M}_{T}$ refers to the empirical measure of
the $\displaystyle N_{min}$–$\displaystyle N_{max}$ particle system at time
$\displaystyle T$, associated to the above branching birth and death process,
and $\displaystyle\mu^{M}_{T}$ refers to the empirical measure of the FV IPS
with $\displaystyle N:=N_{min}=N_{max}$ particles at time $\displaystyle 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 $\displaystyle N_{min}=N_{max}=N$,

$\displaystyle\left\|\frac{m^{M}_{T}(f)}{N}-\nu_{M}(f)\right\|_{2}\leq\frac{(2C_{X}+C)\|f\|_{\infty}}{N^{\alpha/2}}+C_{X}\,e^{-\gamma_{X}T}\,\|f\|_{\infty}.$

for some $\displaystyle\alpha>0$. Using the fact that the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ process is ergodic, we deduce that its
normalised empirical stationary distribution
$\displaystyle\mathcal{X}^{M}_{N}$ (the normalised empirical distribution of a
random vector of particles distributed according to the stationary
distribution of the $\displaystyle N_{min}$–$\displaystyle N_{max}$ process)
satisfies

$\displaystyle\mathcal{X}^{M}_{N}(f)\xrightarrow[N\to+\infty]{L_{2}}\nu_{M}(f),$

for all function $\displaystyle f:\\{1,\ldots,M\\}\to\nu_{M}$. Similarly, the
normalised empirical stationary distribution
$\displaystyle\mathcal{Y}^{M}_{N}$ of the FV particle system satisfies

$\displaystyle\mathcal{Y}^{M}_{N}(f)\xrightarrow[N\to+\infty]{L_{2}}\nu_{M}(f).$

These theoretical convergence results ensure that both algorithms can be used
to approximate the limiting distribution $\displaystyle\nu_{M}$. It remains to
compare their performances via numerical simulations.

##### Numerical comparison of the $\displaystyle N_{min}$–$\displaystyle
N_{max}$ and the FV IPSs.

We now consider the problem of numerically approximating the mean of
$\displaystyle\nu_{M}$, that is $\displaystyle\nu_{M}(f_{M})$ with
$\displaystyle f_{M}(x)=x$ for all $\displaystyle x\in\\{1,\ldots,M\\}$. The
estimator is $\displaystyle\hat{\theta}_{N}=\mathcal{X}^{M}_{N}(f_{M})$ (resp.
$\displaystyle\hat{\theta}_{N}=\mathcal{Y}^{M}_{N}(f_{M})$) for the
$\displaystyle N_{min}$–$\displaystyle N_{max}$ (resp. FV) IPS, and we compare
three important metrics determining the performance of the algorithms:

  * •

the bias of the estimator,
$\displaystyle|\mathbb{E}(\hat{\theta}_{N})-\nu_{M}(f_{M})|$,

  * •

the standard deviation of the estimator,
$\displaystyle\text{Std}(\hat{\theta}_{N})$,

  * •

the number of interaction events, $\displaystyle(A_{T}+B_{T})/T$, which
corresponds to the number $\displaystyle B_{T}$ of selection events in the
$\displaystyle N_{min}$–$\displaystyle N_{max}$ algorithm and the number
$\displaystyle A_{T}$ of resampling events in FV IPS, per unit of time.

The last metric $\displaystyle(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 $\displaystyle N=10$ and then $\displaystyle N=100$, we let the upper
bound of the state space $\displaystyle M$ vary in
$\displaystyle\\{10,100,1000,+\infty\\}$, and present the results in Table 1
for $\displaystyle N=10$, and in Table 2 for $\displaystyle N=100$. For each
metric, a lower number indicates that the algorithm performs better.

$\displaystyle N=10$$\displaystyle\Big{)}$ |  $\displaystyle|\mathbb{E}(\hat{\theta}_{N})-\nu_{M}(f_{M})|$ |  $\displaystyle\text{Std}(\hat{\theta}_{N})$ |  $\displaystyle(A_{T}+B_{T})/T$
---|---|---|---
$\displaystyle M=10$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.08$ |  $\displaystyle 0.30$ |  $\displaystyle 14.0$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.10$ |  $\displaystyle 0.41$ |  $\displaystyle 87.2$
$\displaystyle M=100$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.08$ |  $\displaystyle 0.30$ |  $\displaystyle 14.0$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.20$ |  $\displaystyle 0.51$ |  $\displaystyle 988$
$\displaystyle M=1000$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.08$ |  $\displaystyle 0.30$ |  $\displaystyle 14.0$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.22$ |  $\displaystyle 0.53$ |  $\displaystyle 9989$
$\displaystyle M=+\infty$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.08$ |  $\displaystyle 0.30$ |  $\displaystyle 14.0$
FV |  $\displaystyle\Big{)}$ $\displaystyle*$ |  $\displaystyle*$ |  $\displaystyle*$
Table 1:  For each value of $\displaystyle M\in\\{10,100,1000,+\infty\\}$, we display the bias of the estimator $\displaystyle\hat{\theta}_{N}$, its standard deviation, and the number of events, for the $\displaystyle N_{min}$–$\displaystyle N_{max}$ algorithm (for which $\displaystyle\hat{\theta}_{N}=\mathcal{X}^{M}_{N}(f_{M})$) and for the FV IPS (for which $\displaystyle\hat{\theta}_{N}=\mathcal{Y}^{M}_{N}(f_{M})$) with $\displaystyle N=N_{min}=N_{max}=10$ particles. Note that the FV IPS is not defined when $\displaystyle M=+\infty$. $\displaystyle N=100$$\displaystyle\Big{)}$ |  $\displaystyle|\mathbb{E}(\hat{\theta}_{N})-\nu_{M}(f_{M})|$ |  $\displaystyle\text{Std}(\hat{\theta}_{N})$ |  $\displaystyle(A_{T}+B_{T})/T$
---|---|---|---
$\displaystyle M=10$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.01$ |  $\displaystyle 0.12$ |  $\displaystyle 144$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.02$ |  $\displaystyle 0.18$ |  $\displaystyle 857$
$\displaystyle M=100$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.01$ |  $\displaystyle 0.12$ |  $\displaystyle 144$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.10$ |  $\displaystyle 0.39$ |  $\displaystyle 9866$
$\displaystyle M=1000$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.01$ |  $\displaystyle 0.12$ |  $\displaystyle 144$
FV |  $\displaystyle\Big{)}$ $\displaystyle 0.20$ |  $\displaystyle 0.50$ |  $\displaystyle 99873$
$\displaystyle M=+\infty$ |  $\displaystyle N_{min}$–$\displaystyle N_{max}$ |  $\displaystyle\Big{)}$ $\displaystyle 0.01$ |  $\displaystyle 0.12$ |  $\displaystyle 144$
FV |  $\displaystyle\Big{)}$ $\displaystyle*$ |  $\displaystyle*$ |  $\displaystyle*$
Table 2:  For each value of $\displaystyle M\in\\{10,100,1000,+\infty\\}$, we
display the bias of the estimator $\displaystyle\hat{\theta}_{N}$, its
standard deviation, and the number of events, for the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ algorithm (for which
$\displaystyle\hat{\theta}_{N}=\mathcal{X}^{M}_{N}(f_{M})$) and for the FV IPS
(for which $\displaystyle\hat{\theta}_{N}=\mathcal{Y}^{M}_{N}(f_{M})$) with
$\displaystyle N=N_{min}=N_{max}=100$ particles. Note that the FV IPS is not
defined when $\displaystyle M=+\infty$.

We can see that, for this particular example, the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ algorithm performs better in all situations,
with respect to all metrics. Moreover, the bias of the FV approximation
increases when $\displaystyle M$ increases, while the variance of the
estimator also increases and the computation cost increases. Without surprise,
increasing $\displaystyle M$ does not change the dynamic of the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ algorithm, since very few particles, if any,
reach the boundary $\displaystyle M$ even when $\displaystyle M=10$. Further,
this suggests that it should be possible to extend the $\displaystyle
N_{min}$–$\displaystyle N_{max}$ algorithm to the case of unbounded branching
rates, and this is supported by the numerical simulations when $\displaystyle
M=+\infty$ (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 $\displaystyle N_{min}$–$\displaystyle 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 $\displaystyle\\{0,\dots,n\\}$ for some
$\displaystyle n\geq 1$ with a time inhomogeneous branching rate. More
precisely, let $\displaystyle(X_{t})_{t\geq 0}$ denote a continuous time
random walk on $\displaystyle\mathbb{Z}$ that moves one step to the right with
probability $\displaystyle p\in(0,1)$ and one step to the left otherwise. Each
step is taken after an (independent) exponentially distributed time with
parameter $\displaystyle 1$.

Now let $\displaystyle b,k:\\{0,\dots,n\\}\to[0,\infty)$ be bounded measurable
functions and $\displaystyle\mathcal{E}_{b}$, $\displaystyle\mathcal{E}_{k}$
be exponential random variables, that are independent of each other and
$\displaystyle X$, and define

$\displaystyle\tau_{1}^{b}=\inf\\{t>0:\int_{0}^{t}b(X_{s}){\rm
d}s>\mathcal{E}_{b}\\},$

and

$\displaystyle\tau_{1}^{k}=\inf\\{t>0:\int_{0}^{t}k(X_{s}){\rm
d}s>\mathcal{E}_{k}\\}.$

Further, set $\displaystyle\tau_{1}=\tau_{1}^{b}\wedge\tau_{1}^{k}$. If
$\displaystyle\tau_{1}=\tau_{1}^{b}$, a new particle is added to the system at
position $\displaystyle X_{\tau_{1}^{b}}$, which continues the same behaviour
as the original particle. If $\displaystyle\tau_{1}=\tau_{1}^{k}$, the
particle is removed from the system. We can then define
$\displaystyle\tau_{2}$ to be the first time after $\displaystyle\tau_{1}$
that one of the particles alive at time $\displaystyle\tau_{1}$ branches or is
killed, and this process continues iteratively.

Letting $\displaystyle N_{t}$ denote the number of particles alive at time
$\displaystyle t\geq 0$, and
$\displaystyle\\{X_{t}^{i}\,:\,i=1,\dots,N_{t}\\}$ their positions, the
branching random walk is defined as

$\displaystyle Z_{t}=\sum_{i=1}^{N_{t}}\delta_{X_{t}^{i}},\qquad t\geq 0.$

We now consider a time-inhomogeneous version of the above process. We consider
a sequence of switching times $\displaystyle
0=T_{0}^{\texttt{off}}<T_{1}^{\texttt{on}}<T_{1}^{\texttt{off}}<T_{2}^{\texttt{on}}<T_{2}^{\texttt{off}}<\dots$,
defined via

$\displaystyle
T_{i}^{\texttt{on}}=T_{i-1}^{\texttt{off}}+\mathbf{e}_{i-1}(s_{\texttt{on}}),\qquad
i\geq 1,$

and

$\displaystyle
T_{i}^{\texttt{off}}=T_{i}^{\texttt{on}}+\mathbf{e}_{i}(s_{\texttt{off}}),\qquad
i\geq 1,$

where $\displaystyle\mathbf{e}_{i}(s_{\texttt{off}})$ and
$\displaystyle\mathbf{e}_{i}(s_{\texttt{on}})$ are independent (of each other
and everything else) exponential random variables with rate $\displaystyle
s_{\texttt{off}}$ and $\displaystyle s_{\texttt{on}}$, respectively. For
$\displaystyle T_{i-1}^{\texttt{off}}\leq t<T_{i}^{\texttt{on}}$,
$\displaystyle Z_{t}$ evolves as above but with $\displaystyle b=0$, i.e. each
particle moves as a random walk until it is (possibly) killed, until time
$\displaystyle T_{i}^{\texttt{on}}$. At time $\displaystyle
T_{i}^{\texttt{on}}$, we sample from an exponential random variable with rate
$\displaystyle B\in(0,\infty)$ and set $\displaystyle b$ to be this value. For
$\displaystyle T_{i}^{\texttt{on}}\leq t<T_{i}^{\texttt{off}}$, $\displaystyle
Z_{t}$ then evolves as the branching random walk with branching rate
$\displaystyle b$ for each particle. At time $\displaystyle
T_{i}^{\texttt{off}}$, we again set the branching rate to be $\displaystyle
0$, and continue this process until some terminal time $\displaystyle 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 $\displaystyle N_{min}$–$\displaystyle N_{max}$ process due to
large deviations effects. We refer the reader to the next section for a
discussion and possible solution to this.

####  3.3.2 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 $\displaystyle Q$ for large times. A natural conjecture
is that for large times $\displaystyle T$, $\displaystyle
m_{0}Q_{T}f\approx\langle
m_{0},\psi\rangle\langle\varphi,f\rangle\mathrm{e}^{\lambda T}$, where
$\displaystyle\varphi$ and $\displaystyle\psi$ are the left eigenmeasure and
right-eigenfunction respectively of the operator $\displaystyle Q$, and
$\displaystyle\lambda$ is the Lyapunov exponent for $\displaystyle Q$. One can
try to estimate $\displaystyle\lambda$ by computing the quantity
$\displaystyle m_{0}Q_{T}f$ for large times $\displaystyle T$.

By Theorem 6, one method for estimating $\displaystyle\lambda$ would be
through equation (6), so that we can simulate the
$\displaystyle(m_{T},\Pi_{T}^{A},\Pi_{T}^{B})$-system and estimate
$\displaystyle\lambda$ by:

$\hat{\lambda}=\frac{1}{T}\log\left(\Pi_{T}^{A}\Pi_{T}^{B}m_{T}(\mathbf{1})\right)=\frac{1}{T}\log\left(\Pi_{T}^{A}\Pi_{T}^{B}\right),$
(33)

for $\displaystyle T$ large. In practice, however, this is numerically a poor
choice of estimator, due to large-deviation effects. The presence of the
‘exponential’ functionals, $\displaystyle\Pi_{T}^{A},\Pi_{T}^{B}$ mean that
significant contributions to the expectation will be made by paths of the
$\displaystyle N_{min}$–$\displaystyle N_{max}$ process which have probability
$\displaystyle~{}e^{-\delta T}$, and so infeasibly many simulations will be
needed to compute an unbiased estimate for $\displaystyle\hat{\lambda}$. A
simple solution is to replace the estimator by

$\bar{\lambda}=\frac{1}{\delta
t}\log\left[\frac{1}{n}\sum_{i<n}\left(\frac{\Pi_{t_{i+1}}^{A}\Pi_{t_{i+1}}^{B}}{\Pi_{t_{i}}^{A}\Pi_{t_{i}}^{B}}\right)\right],$
(34)

where $\displaystyle t_{i}=T\frac{i}{N}$ and $\displaystyle\delta
t=\frac{T}{N}$. In many examples, under conditions which ensure that the
process is ergodic and thus converges to a stationary distribution, as
$\displaystyle T\to\infty$, for fixed $\displaystyle\delta t$, the estimator
$\displaystyle\bar{\lambda}$ can be shown to converge to
$\displaystyle\lambda$ (e.g. [24]).

In some cases the expression in (34) may not be valid and we can only use
estimates of the form in (33). Examples may be when the process is not
ergodic, or the system is not time-homogenous, such as the example given in
Remark 22. In such circumstances, we propose a novel Monte Carlo method to
compute estimates of $\displaystyle\lambda$. Specifically, we note that the
$\displaystyle N_{min}$–$\displaystyle N_{max}$ process
$\displaystyle\mathcal{X}_{t}:=\\{X_{t}^{1},\ldots,X_{t}^{N_{t}}\\}$ defines a
continuous time Markov process, while the processes $\displaystyle A_{t}$ and
$\displaystyle B_{t}$ corresponding to the number of selection and resampling
events, are themselves additive functionals of the Markov process. Estimation
of $\displaystyle\lambda$ 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, [21, 20, 19, 18]. This suggests an algorithm as follows,
assuming that $\displaystyle b^{i}$, $\displaystyle\kappa^{i}$, $\displaystyle
p^{i}$ and $\displaystyle q^{i}$ have been chosen:

  1. 1.

Fix $\displaystyle T,N^{*},\delta t,N_{0},x_{0}$. Initiate with $\displaystyle
t=0$, and $\displaystyle\mathcal{X}^{1},\ldots,\mathcal{X}^{N^{*}}$ each
containing $\displaystyle N_{0}$ particles of initial value $\displaystyle x$.
Set $\displaystyle W=1$.

  2. 2.

Run each $\displaystyle\mathcal{X}^{i}$ from time $\displaystyle t$ to time
$\displaystyle t^{\prime}=\min\\{t+\delta t,T\\}$ according to the
$\displaystyle N_{min}-N_{max}$ algorithm above, setting $\displaystyle
w_{i}:=\frac{\Pi_{t^{\prime}}^{A^{i}}\Pi_{t^{\prime}}^{B^{i}}}{\Pi_{t}^{A^{i}}\Pi_{t}^{B^{i}}}$.

  3. 3.

Set $\displaystyle
W:=W\times\left(\frac{1}{N*}\sum_{i=1}^{N^{*}}w_{i}\right)$.

  4. 4.

If $\displaystyle t^{\prime}<T$, set $\displaystyle t=t^{\prime}$, resample
each $\displaystyle\mathcal{X}^{i}$ from the existing particles
$\displaystyle\\{\mathcal{X}^{1},\ldots,\mathcal{X}^{N^{*}}\\}$ with
probabilities proportional to $\displaystyle w_{1},\ldots,w_{N^{*}}$,
independently (multinomial resampling). Return to step 2.

  5. 5.

If $\displaystyle t^{\prime}=T$, stop, and return

$\hat{\lambda}:=\frac{1}{T}\log\left(W\right).$

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
$\displaystyle\lambda$.

(a) Single Trajectory Simulation

(b) Particle Filter Simulation

Figure 2: We use two methods described above to estimate the eigenvalue
$\displaystyle\lambda$. The first simulation (left) estimates
$\displaystyle\lambda$ using (33), based on a single trajectory run up to time
$\displaystyle T=4000$. The true value of $\displaystyle\lambda$ 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), $\displaystyle\lambda$
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
$\displaystyle 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.

##  4 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 $\displaystyle\mathbb{M}$ and
$\displaystyle\mathcal{M}$ such that

$\displaystyle\displaystyle\nu_{t}^{f}-\nu_{0}^{f}$
$\displaystyle\displaystyle=\Pi_{t}^{A}\Pi_{t}^{B}(\mathbb{M}_{t}-\mathbb{M}_{\tau_{C_{t}}})+\sum_{n=1}^{A_{t}}\Pi_{\rho_{n}}^{A}\Pi_{\rho_{n}}^{B}\left(\mathcal{M}_{\rho_{n}}-\mathcal{M}_{\rho_{n}-}\right)+\sum_{n=1}^{B_{t}}\Pi_{\sigma_{n}}^{A}\Pi_{\sigma_{n}}^{B}(\mathcal{M}_{\sigma_{n}}-\mathcal{M}_{\sigma_{n}-})$
$\displaystyle\displaystyle\phantom{\Pi_{t}^{A}\Pi_{t}^{B}(\mu_{t}-\mu_{\tau_{C_{t}}}}\phantom{+\sum_{n=1}^{B_{t}}\Pi^{A}_{\sigma_{n-1}}\Pi^{B}_{\rho_{n-1}}(\mathbb{M}_{\rho_{n}}-\mathbb{M}_{\rho_{n-1}})}+\sum_{n=1}^{C_{t}}\Pi^{A}_{\tau_{n-1}}\Pi^{B}_{\tau_{n-1}}(\mathbb{M}_{\tau_{n}}-\mathbb{M}_{\tau_{n-1}}),$
(35)

where $\displaystyle\sigma_{n}$ and $\displaystyle\rho_{n}$ denote the
$\displaystyle n^{th}$ selection and resampling events, respectively, and

$\displaystyle\nu_{t}^{f}=\Pi_{t}^{A}\Pi_{t}^{B}\sum_{i\in\bar{S}_{t}}Q_{T-t}f(X_{t}^{i}).$

The martingale $\displaystyle\mathbb{M}$ compensates the evolution of the
branching process between resampling/selection times. The martingale
$\displaystyle\mathcal{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 $\displaystyle L^{2}$ norms of the
martingale increments in order to obtain (7).

######  Lemma 23.

For all $\displaystyle\ell\geq 1$, for all $\displaystyle T\geq 0$,

$\displaystyle\mathbb{E}\left(\left(\Pi_{T}^{B}\right)^{\ell}\right)\leq\exp(c_{\ell}\|b\|_{\infty}T),$

where

$\displaystyle c_{\ell}:=\sup_{x\geq 2}x\left((1+1/x)^{\ell}-1\right).$

In addition, denoting by $\displaystyle\beta_{T}$ the total number of
branching events (with or without selection) up to time $\displaystyle T$, we
have

$\displaystyle\mathbb{E}\left(\beta_{T}\,\left(\Pi_{T}^{B}\right)^{2}\right)\leq
N_{0}\exp((3+c_{4})\|b\|_{\infty}T/2).$

The above estimates are the main reason why we require $\displaystyle 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 $\displaystyle\mathbb{M}$. For all $\displaystyle
i\in\mathbb{N}$, we define respectively the birth time
$\displaystyle\texttt{B}_{i}$ and the death time $\displaystyle\texttt{D}_{i}$
of the particle with index $\displaystyle i$ as

$\displaystyle\displaystyle\texttt{B}_{i}=\inf\\{t\geq 0,\
i\in\bar{S}_{t}\\}\text{ and }\texttt{D}_{i}=\inf\\{t\geq\texttt{B}_{i},\
i\notin\bar{S}_{t}\\}.$

Fix $\displaystyle T>0$ and $\displaystyle t\in[0,T]$. For all particles
$\displaystyle i\in\mathbb{N}$, we set

$\mathbb{M}^{i,n}_{t}=\begin{cases}0,\quad\text{if
}t\leq\texttt{B}_{i}\vee\tau_{n},\\\
Q_{T-t}f(X^{i}_{t})-Q_{T-\texttt{B}_{i}\vee\tau_{n}}f(X_{\texttt{B}_{i}\vee\tau_{n}}^{i}),\quad\text{if
}t\in[\texttt{B}_{i}\vee\tau_{n},\texttt{D}_{i}\wedge\tau_{n+1}),\\\
(1-\mathbf{1}_{D^{n}(i)}+\mathbf{1}_{B^{n}(i)})Q_{T-\texttt{D}_{i}\wedge\tau_{n+1}}f(X_{\texttt{D}_{i}\wedge\tau_{n+1}}^{i})-Q_{T-\texttt{B}_{i}\vee\tau_{n}}f(X^{i}_{\texttt{B}_{i}\vee\tau_{n}}),\quad\text{if
}t\geq\texttt{D}_{i}\wedge\tau_{n+1},\end{cases}$ (36)

where $\displaystyle B^{n}(i)$ denotes the event “$\displaystyle\tau_{n+1}$ is
a branching time for the $\displaystyle i$-th particle” and $\displaystyle
D^{n}(i)$ denotes the event “$\displaystyle\tau_{n+1}=\texttt{D}_{i}$ is a
soft killing time for the $\displaystyle i$-th particle”.

By definition of the process between the time intervals
$\displaystyle[\tau_{n},\tau_{n+1})$, to show that $\displaystyle\mathbb{M}$
is indeed a martingale, it is sufficient to prove that

$\displaystyle\mathbb{M}^{1,0}_{t}:=\begin{cases}Q_{T-t}f(X^{1}_{t})-Q_{T}f(X^{1}_{0}),\quad\text{if
}t\in[0,\tau_{1}),\\\
(1-\mathbf{1}_{D(1)}+\mathbf{1}_{B(1)})Q_{T-\tau_{1}}f(X_{\tau_{1}})-Q_{T}f(X_{0}),\quad\text{if
}t\geq\tau_{1},\end{cases}$

is a martingale. We observe that, by construction,
$\displaystyle\tau_{1}=\tau^{\partial}\wedge\tau^{\kappa}\wedge\tau^{b}$ where

$\displaystyle\displaystyle\tau^{b}$
$\displaystyle\displaystyle=\inf_{j\in\bar{S}_{0}}\inf\\{t\geq 0,\
\int_{0}^{t}b^{j}((X^{i}_{u})_{i\in\bar{S}_{u}})\,\mathrm{d}u\geq e^{b}\\},$
$\displaystyle\displaystyle\tau^{\kappa}$
$\displaystyle\displaystyle=\inf_{j\in\bar{S}_{0}}\inf\\{t\geq 0,\
\int_{0}^{t}\kappa^{j}((X^{i}_{u})_{i\in\bar{S}_{u}})\,\mathrm{d}s\geq
e^{\kappa}\\},$ $\displaystyle\displaystyle\tau^{\partial}$
$\displaystyle\displaystyle=\inf_{j\in\bar{S}_{0}}\inf\\{t\geq 0\
X^{j}_{t}\in\partial\\}.$

Then, for all $\displaystyle s,t\geq 0$ such that $\displaystyle s+t\leq T$,
setting
$\displaystyle\mathbb{E}^{t}=\mathbb{E}[\cdot\mid\sigma(\bar{S}_{u},\,\bar{X}_{u},0\leq
u\leq t)]$ and using the fact that $\displaystyle Q_{u}f\equiv 0$ on
$\displaystyle\partial$ for all $\displaystyle u\geq 0$, we have

$\displaystyle\displaystyle\mathbb{E}^{t}(\mathbb{M}^{1,0}_{t+s})=\mathbb{E}^{t}$
$\displaystyle\displaystyle\left[\left(Q_{T-(t+s)}f(X^{1}_{t+s})-Q_{T}f(X_{0})\right)\mathbf{1}_{t+s<\tau_{1}}\right]$
$\displaystyle\displaystyle\qquad+\mathbb{E}_{x}^{t}\left[\left((1-\mathbf{1}_{D(1)}+\mathbf{1}_{B(1)})Q_{T-\tau_{1}}f(X_{\tau_{1}})-Q_{T}f(X_{0})\right){\mathbf{1}}_{t<\tau_{1}\leq
t+s}\right]$
$\displaystyle\displaystyle\qquad+\left((1-\mathbf{1}_{D(1)}+\mathbf{1}_{B(1)})Q_{T-\tau_{1}}f(X_{\tau_{1}})-Q_{T}f(X_{0})\right){\mathbf{1}}_{\tau_{1}\leq
t}$ (37)

Using the fact that the particles are independent copies of $\displaystyle X$
before time $\displaystyle\tau_{1}$, and the Markov property at time
$\displaystyle t$, we obtain

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}\left[Q_{T-(t+s)}f(X^{1}_{t+s})\mathbf{1}_{t+s<\tau_{1}}\right]$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[Q_{T-(t+s)}f(Y^{1}_{t+s})e^{-\int_{t}^{t+s}h_{u}\,{\rm
d}u}{\mathbf{1}}_{t+s<\tau_{\partial}}\right]\mathbf{1}_{t<\tau_{1}},$

where $\displaystyle\texttt{E}^{t}$ is the expectation on a probability space
$\displaystyle\Omega^{\prime}$ such that $\displaystyle(Y^{i}_{u})_{u\geq t}$,
$\displaystyle i\in\bar{S}_{0}$, are independent copies of $\displaystyle X$,
starting from $\displaystyle X^{i}_{t}$ at time $\displaystyle t$, and where

$\displaystyle
h_{u}:=\sum_{j\in\bar{S}_{0}}b^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}})+\kappa^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}}).$

Then, using the Markov property at time $\displaystyle t+s$, we obtain

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}\left[Q_{T-(t+s)}f(X^{1}_{t+s})\mathbf{1}_{t+s<\tau_{1}}\right]$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[e^{\int_{t+s}^{T}b(Y^{1}_{u})-\kappa(Y^{1}_{u})\,{\rm
d}u}f(Y^{1}_{T})e^{-\int_{t}^{t+s}h_{u}\,{\rm
d}u}{\mathbf{1}}_{t+s<\tau_{\partial}}\right]\mathbf{1}_{t<\tau_{1}}$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[e^{\int_{t}^{T}b(Y^{1}_{u})-\kappa(Y^{1}_{u})\,{\rm
d}u}f(Y^{1}_{T})e^{-\int_{t}^{t+s}h^{1}_{u}\,{\rm
d}u}{\mathbf{1}}_{t+s<\tau_{\partial}}\right]\mathbf{1}_{t<\tau_{1}},$

where $\displaystyle
h^{1}_{u}:=2b^{1}((Y^{i}_{u})_{i\in\bar{S}_{0}})+\sum_{j\in\bar{S}_{0},j\neq
1}b^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}})+\kappa^{j}((Y^{i}_{u})_{i\in\bar{S}_{0}})$.
Similarly, since $\displaystyle\kappa^{1}$ is the intensity of the soft
killing time for the first particle, and $\displaystyle b^{1}$ is the
intensity of its branching time, we obtain

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}\left[(1-{\mathbf{1}}_{D(1)}+{\mathbf{1}}_{B(1)})Q_{T-\tau_{1}}f(X^{1}_{\tau_{1}})\mathbf{1}_{t<\tau_{1}\leq
t+s}\right]$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[\int_{t}^{(t+s)\wedge\tau_{\partial}}(h_{u}-\kappa^{1}((Y^{i}_{u})_{i\in\bar{S}_{0}})+b^{1}((Y^{i}_{u})_{i\in\bar{S}_{0}}))Q_{T-u}f(Y^{1}_{u})e^{-\int_{t}^{u}h_{v}\,\mathrm{d}v}\mathrm{d}u\right]\,{\mathbf{1}}_{t<\tau_{1}}$
$\displaystyle\displaystyle\qquad+\texttt{E}^{t}\left[(1-{\mathbf{1}}_{Y^{1}_{\tau_{\partial}}\in{\partial}})e^{-\int_{t}^{\tau_{\partial}}h_{u}\,\mathrm{d}u}{\mathbf{1}}_{\tau_{\partial}\leq
t+s}Q_{T-\tau_{\partial}}f(Y^{1}_{\tau_{\partial}})\right]\,{\mathbf{1}}_{t<\tau_{1}}.$

Hence, using the fact that $\displaystyle Y$ satisfies the strong Markov
property at time $\displaystyle\tau_{\partial}$, and the fact that
$\displaystyle Q_{T-\tau_{\partial}}f(Y^{1}_{\tau_{\partial}})=0$ on the event
$\displaystyle\\{Y^{1}_{\tau_{\partial}}\in{\partial}\\}$, we deduce that

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}\left[(1-{\mathbf{1}}_{D(1)}+{\mathbf{1}}_{B(1)})Q_{T-\tau_{1}}f(X^{1}_{\tau_{1}})\mathbf{1}_{t<\tau_{1}\leq
t+s}\right]$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[\int_{t}^{(t+s)\wedge\tau_{\partial}}h^{1}_{u}\,e^{\int_{u}^{T}(b(Y^{1}_{v})-\kappa(Y^{1}_{v}))\,\mathrm{d}v}\,f(Y^{1}_{T})e^{-\int_{t}^{u}h_{v}\,\mathrm{d}v}\,\mathrm{d}u\right]\,{\mathbf{1}}_{t<\tau_{1}}$
$\displaystyle\displaystyle\qquad+\texttt{E}^{t}\left[e^{-\int_{t}^{\tau_{\partial}}h_{u}\,\mathrm{d}u}{\mathbf{1}}_{\tau_{\partial}\leq
t+s}e^{\int_{\tau_{\partial}}^{T}(b(Y^{1}_{v})-\kappa(Y^{1}_{v}))\,\mathrm{d}v}f(Y^{1}_{T})\right]\,{\mathbf{1}}_{t<\tau_{1}}$
$\displaystyle\displaystyle=\texttt{E}^{t}\left[e^{\int_{t}^{T}(b(Y^{1}_{v})-\kappa(Y^{1}_{v}))\,\mathrm{d}v}\,f(Y^{1}_{T})\,\left(\int_{t}^{(t+s)\wedge\tau_{\partial}}h^{1}_{u}\,e^{-\int_{t}^{u}h^{1}_{v}\,\mathrm{d}v}\mathrm{d}u+e^{-\int_{t}^{\tau_{\partial}}h^{1}_{v}\,\mathrm{d}v}{\mathbf{1}}_{\tau_{\partial}\leq
t+s}\right)\right]\,{\mathbf{1}}_{t<\tau_{1}},$

where we used the assumption $\displaystyle
b(Y^{1}_{v})-\kappa(Y^{1}_{v})=b^{1}((Y^{i}_{u})_{i\in\bar{S}_{0}})-\kappa^{1}((Y^{i}_{u})_{i\in\bar{S}_{0}})$.
But

$\displaystyle e^{-\int_{t}^{t+s}h^{1}_{u}\,{\rm
d}u}{\mathbf{1}}_{t+s<\tau_{\partial}}+\int_{t}^{(t+s)\wedge\tau_{\partial}}h^{1}_{u}\,e^{-\int_{t}^{u}h^{1}_{v}\,\mathrm{d}v}\mathrm{d}u+e^{-\int_{t}^{\tau_{\partial}}h^{1}_{v}\,\mathrm{d}v}{\mathbf{1}}_{\tau_{\partial}\leq
t+s}=1$

Hence we deduce from (4) that

$\displaystyle\displaystyle\mathbb{E}_{x}^{t}(\mathbb{M}^{1,0}_{t+s})$
$\displaystyle\displaystyle=(1-\mathbf{1}_{D(1)}+\mathbf{1}_{B(1)})\left(Q_{T-\tau_{1}}f(X_{\tau_{1}})-Q_{T}f(X_{0})\right){\mathbf{1}}_{\tau_{1}\leq
t}$
$\displaystyle\displaystyle\qquad+\left(\texttt{E}^{t}\left[e^{\int_{t}^{T}(b(Y^{1}_{v})-\kappa(Y^{1}_{v}))\,\mathrm{d}v}\,f(Y^{1}_{T})\right]-Q_{T}f(x)\right){\mathbf{1}}_{t<\tau_{1}}$
$\displaystyle\displaystyle=(1-\mathbf{1}_{D(1)}+\mathbf{1}_{B(1)})Q_{T-\tau_{1}}f(X_{\tau_{1}}){\mathbf{1}}_{\tau_{1}\leq
t}+Q_{T-t}f(X^{1}_{t})\,{\mathbf{1}}_{t<\tau_{1}}-Q_{T}f(x)$
$\displaystyle\displaystyle=\mathbb{M}^{1,0}_{t}.$

This concludes Step 1.

Step 2. A first martingale decomposition.

Fix $\displaystyle t\geq 0$ and recall that $\displaystyle C_{t}=\inf\\{n\geq
0:t<\tau_{n+1}\\}$. In what follows, we denote by
$\displaystyle\rho_{1},\rho_{2},\ldots$ the sequence of times at which a
resampling event occurs, and by $\displaystyle\sigma_{1},\sigma_{2},\ldots$
the sequence of times at which a selection event occurs. Note that, if
$\displaystyle\tau_{n}\in\\{\rho_{1},\rho_{2},\ldots\\}$ or
$\displaystyle\tau_{n}\in\\{\sigma_{1},\sigma_{2},\ldots\\}$, then
$\displaystyle N_{n}=N_{n-1}\in\\{N_{min},\ldots,N_{max}\\}$.

We will show that, for all $\displaystyle t\geq 0$, we have

$\displaystyle\displaystyle\mathbb{M}_{t}:=\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}=$
$\displaystyle\displaystyle\sum_{i\in\bar{S}_{t}}Q_{T-t}f(X_{t}^{i})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X_{0}^{i})$
$\displaystyle\displaystyle\qquad-\sum_{n=1}^{A_{t}}Q_{T-\rho_{n}}f(X_{\rho_{n}}^{i_{n}})+\sum_{n=1}^{B_{t}}Q_{T-\sigma_{n}}f(X_{\sigma_{n}}^{j_{n-1}}),$
(38)

where we recall that $\displaystyle A_{t}$ is the number of resampling events
up to time $\displaystyle t$, $\displaystyle B_{t}$ is the number of selection
events up to time $\displaystyle t$, and where we have set $\displaystyle
i_{n}$ to denote the index of the particle added at time
$\displaystyle\rho_{n}$ and $\displaystyle j_{n-1}$ is the particle killed
(due to a selection event) at time $\displaystyle\sigma_{n}$ (as above,
$\displaystyle X_{\rho_{n}}^{j_{n-1}}$ denotes the position of the particle
before the selection event).

We will show that (38) holds for $\displaystyle n=0,1$, since the general
formula then follows using similar arguments.

Noting that for all $\displaystyle t\in[\tau_{0},\tau_{1})$, $\displaystyle
C_{t}=A_{t}=B_{t}=0$, in this case we have

$\displaystyle\displaystyle\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}=\sum_{i\in\bar{S}_{0}}\mathbb{M}_{t}^{i,0}=\sum_{i\in\bar{S}_{0}}Q_{T-t}f(X_{t}^{i})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X_{0}^{i}),$

as required.

Let us now describe the evolution of the martingale at time
$\displaystyle\tau_{1}$, depending on the type of event that occurs.

Case $\displaystyle\tau_{1}=\rho_{1}$. In this case, for all $\displaystyle
t\in[\tau_{1},\tau_{2})$, we have $\displaystyle A_{t}=1$ and $\displaystyle
B_{t}=0$, and
$\displaystyle\bar{S}_{t}=\\{1+\max\bar{S}_{0}\\}\cup\bar{S}_{0}\setminus\\{i^{\prime}_{0}\\}$,
where $\displaystyle i^{\prime}_{0}$ is the particle removed at time
$\displaystyle\tau_{1}$. Thus

$\displaystyle\displaystyle\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{0}}\mathbb{M}_{t}^{i,0}+\sum_{i\in\bar{S}_{\tau_{1}}}\mathbb{M}_{t}^{i,1}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{0}}\left[(1-{\mathbf{1}}_{i=i^{\prime}_{0}})Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i})-Q_{T}f(X_{0}^{i})\right]+\sum_{i\in\bar{S}_{\tau_{1}}}\big{[}Q_{T-t}f(X^{i}_{t})-Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i})\big{]}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{\tau_{1}}}Q_{T-t}f(X_{t}^{i})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X_{0}^{i})-Q_{T-\tau_{1}}f(X_{\tau_{1}}^{1+\max\bar{S}_{0}})$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{t}}Q_{T-t}f(X_{t}^{i})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X_{0}^{i})-Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i_{1}}),$

as claimed.

Case $\displaystyle\tau_{1}=\sigma_{1}$. Similarly, if
$\displaystyle\tau_{1}=\sigma_{1}$, then, for all $\displaystyle
t\in[\tau_{1},\tau_{2})$, we have $\displaystyle A_{t}=0$ and $\displaystyle
B_{t}=1$, and
$\displaystyle\bar{S}_{t}=\\{1+\max\bar{S}_{0}\\}\cup\bar{S}_{0}\setminus\\{j_{0}\\}$.
Denoting by $\displaystyle j^{\prime}_{0}$ the particle duplicated at time
$\displaystyle\tau_{1}$, we have

$\displaystyle\displaystyle\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{0}}\mathbb{M}_{t}^{i,0}+\sum_{i\in\bar{S}_{\tau_{1}}}\mathbb{M}_{t}^{i,1}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{0}}\left[(1+{\mathbf{1}}_{i=j^{\prime}_{0}})Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i})-Q_{T}f(X_{0}^{i})\right]+\sum_{i\in\bar{S}_{\tau_{1}}}\big{[}Q_{T-t}f(X^{i}_{t})-Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i})\big{]}$
$\displaystyle\displaystyle=\sum_{i\in S_{1}}Q_{T-t}f(X^{i}_{t})-\sum_{i\in
S_{0}}Q_{T}f(X_{0}^{i})+Q_{T-\tau_{1}}f(X^{j_{0}}_{\tau_{1}})$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{t}}Q_{T-t}f(X^{i}_{t})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X_{0}^{i})+Q_{T-\tau_{1}}f(X^{j_{0}}_{\tau_{1}}).$

Case $\displaystyle\tau_{1}\notin\\{\rho_{1},\sigma_{1}\\}$. If
$\displaystyle\tau_{1}\notin\\{\rho_{1},\sigma_{1}\\}$, we have, for all
$\displaystyle t\in[\tau_{1},\tau_{2})$, $\displaystyle A_{t}=B_{t}=0$. If
$\displaystyle\tau_{1}$ is a branching event of the particle with index
$\displaystyle j^{\prime}_{0}\in S_{0}$, then
$\displaystyle\bar{S}_{t}=\\{1+\max\bar{S}_{0}\\}\cup\bar{S}_{0}$

$\displaystyle\displaystyle\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{0}}\left[(1+{\mathbf{1}}_{i=j^{\prime}_{0}})Q_{T-\tau_{1}}f(X^{i}_{\tau_{1}})-Q_{T}f(X_{0}^{i})\right]+\sum_{i\in\bar{S}_{\tau_{1}}}\left[Q_{T-t}f(X^{i}_{t})-Q_{T-\tau_{1}}f(X_{\tau_{1}}^{i})\right]$
$\displaystyle\displaystyle=\sum_{i\in\bar{S}_{\tau_{1}}}Q_{T-t}f(X^{i}_{t})-\sum_{i\in\bar{S}_{0}}Q_{T}f(X^{i}_{0}),$

where we used the fact that $\displaystyle
X^{j^{\prime}_{0}}_{\tau_{1}}=X^{1+\max\bar{S}_{0}}_{\tau_{1}}$.

Continuing in this manner, we obtain (38).

Step 3. A second martingale decomposition.

For all time $\displaystyle t\geq 0$ and all bounded measurable function
$\displaystyle f:E\to\mathbb{R}$, we denote the integral of $\displaystyle f$
with respect to the occupation measure of the particle system at time
$\displaystyle t$ by

$\mu_{t}^{f}\coloneqq\sum_{i\in\bar{S}_{t}}Q_{T-t}f(X_{t}^{i}).$ (39)

Further, recall the definitions of $\displaystyle\Pi_{t}^{A}$,
$\displaystyle\Pi_{t}^{B}$ and
$\displaystyle\nu_{t}^{f}=\Pi_{t}^{A}\Pi_{t}^{B}\mu_{t}^{f}$.

In this step, we prove that equation (35) holds with $\displaystyle\mathbb{M}$
defined as in the previous steps, and $\displaystyle\mathcal{M}$ to be defined
shortly. To this end, we first note that from Step 2 we obtain

$\displaystyle\displaystyle\mu_{t}^{f}-\mu_{0}^{f}=\sum_{n=0}^{C_{t}}\sum_{i\in\bar{S}_{t}}\mathbb{M}_{t}^{i,n}$
$\displaystyle\displaystyle+\sum_{n=1}^{A_{t}}Q_{T-\rho_{n}}f(X_{\rho_{n}}^{i_{n-1}})-\sum_{n=1}^{B_{t}}Q_{T-\sigma_{n}}f(X_{\sigma_{n}}^{j_{n-1}}).$
(40)

For $\displaystyle n\in\mathbb{N}$, consider the following martingale
increments. We define

$\mathcal{M}_{\rho_{n}}-\mathcal{M}_{\rho_{n}-}\coloneqq
Q_{T-\rho_{n}}f(X_{\rho_{n}}^{i_{n-1}})-\frac{1}{\bar{N}_{\rho_{n}}-1}\sum_{i\in\bar{S}_{\rho_{n}-}\setminus\\{i^{\prime}_{n-1}\\}}Q_{T-\rho_{n}}f(X_{\rho_{n}}^{i}),$
(41)

where $\displaystyle i_{n-1}$ is the index of the particle duplicated at time
$\displaystyle\rho_{n}$ and $\displaystyle i^{\prime}_{n-1}$ is the index of
the particle removed at time $\displaystyle\rho_{n}$ (note also that
$\displaystyle\bar{N}_{\rho_{n}}=\bar{N}_{\rho_{n}-}=|\bar{S}_{\rho_{n}-}|$).
We also define

$\mathcal{M}_{\sigma_{n}}-\mathcal{M}_{\sigma_{n}-}\coloneqq\frac{1}{\bar{N}_{\sigma_{n}}+1}\left(\sum_{i\in\bar{S}_{\sigma_{n}-}}Q_{T-\sigma_{n}}f(X_{\sigma_{n}}^{i})+Q_{T-\sigma_{n}}f(X_{\sigma_{n}}^{j^{\prime}_{n-1}})\right)-Q_{T-\sigma_{n}}f(X_{\sigma_{n}}^{j_{n-1}}),$
(42)