11footnotetext: Madynes team, INRIA Nancy – Grand Est, IECL – UMR 7503, Nancy-
Université, Campus scientifique, B.P. 70239, 54506 Vandœuvre-lès-Nancy Cedex,
France22footnotetext: LORIA – UMR 7503, Nancy-Université, Campus scientifique,
B.P. 70239, 54506 Vandœuvre-lès-Nancy Cedex, France, E-mail:
<EMAIL_ADDRESS>
# Approximations of the allelic frequency spectrum in general supercritical
branching populations

Benoit Henry1,2

###### Abstract

We consider a general branching population where the lifetimes of individuals
are i.i.d. with arbitrary distribution and where each individual gives birth
to new individuals at Poisson times independently from each other. In
addition, we suppose that individuals experience mutations at Poissonian rate
$\theta$ under the infinitely many alleles assumption assuming that types are
transmitted from parents to offspring. This mechanism leads to a partition of
the population by type, called the allelic partition. The main object of this
work is the frequency spectrum $A(k,t)$ which counts the number of families of
size $k$ in the population at time $t$. The process $(A(k,t),\
t\in\mathbb{R}_{+})$ is an example of non-Markovian branching process
belonging to the class of general branching processes counted by random
characteristics. In this work, we propose methods of approximation to replace
the frequency spectrum by simpler quantities. Our main goal is study the
asymptotic error made during these approximations through central limit
theorems. In a last section, we perform several numerical analysis using this
model, in particular to analyze the behavior of one of these approximations
with respect to Sabeti’s Extended Haplotype Homozygosity [18].

MSC 2000 subject classifications: Primary 60J80; secondary 92D10, 60J85,
60G51, 60K15, 60F05.

Key words and phrases. branching process – splitting tree – Crump–Mode–Jagers
process – linear birth–death process – Central Limit Theorem.

##  1 Introduction

In this paper, we consider a general branching population where the lifetimes
of the individuals and their reproductions processes are i.i.d. Moreover, we
assume that their lifetimes are distributed according to an arbitrary
probability distribution $\mathbb{P}_{V}$ and that the births occur, during
their lifetime, according to a Poisson process with rate $b$. The tree
underlying this dynamics is called a splitting tree. This class of random
trees was introduced in [11] by Geiger and Kersting and has been widely
studied in the last decade [14, 15, 16].

We suppose, in addition, that mutations occur on individuals and that each new
mutation confers to its holder a brand new type (i.e. never seen in the
population): this is the _infinitely many alleles_ assumption. This allows
modeling the occurrence of a new type in a population (such as a new species
or a new phenotype in a given species). We also suppose that every individual
inherits the type of its parent. This model leads to a partition of the
population by types. The frequency spectrum of the population alive at time
$t$ is defined as the sequence of number $\left(A(k,t)\right)_{k\geq 1}$
where, for each $k$, $A(k,t)$ is the number of families of size $k$ in the
population. The famous example of Ewens sampling formula gives explicit
expression for the law of the frequency spectrum [9] when the genealogy is
given by the Kingman’s coalescent. Other works studied similar quantities in
the case of Galton-Waston branching processes (see [4] or [12]). In our model,
the frequency spectrum has also been widely studied in the past [6, 7, 8, 5].

Another object of interest is the process $(N_{t},\ t\in\mathbb{R}_{+})$ which
counts the number of living individuals in the population at a given time $t$.
This process is known as binary homogeneous Crump-Mode-Jagers process. One of
the main result of the theory of such process is the law of large number which
gives in our particular case that $e^{-\alpha t}N_{t}$ converges almost surely
to a random variable $\mathcal{E}$ which is exponential conditionally on non-
extinction (for some positive constant $\alpha$).

As for $e^{-\alpha t}N_{t}$, it is also known that the quantities $e^{-\alpha
t}A(k,t)$ converge almost surely to $c_{k}\mathcal{E}$, where $c_{k}$ is an
explicit constant. This result can be easily obtained by conjunction of the
works of [6] and [17] using the theory of general branching processes counted
by random characteristics (a complete statement can be found in [8]). An
alternative proof avoiding the use of the general branching processes theory
can be found in [5].

It appears that the frequency spectrum $(A(k,t))_{k\geq 1}$ is a quantity
which is hard to manipulate from the probabilistic point of view (see [6, 7,
5]). This implies that such a model is inconvenient for practical
applications. In this work we propose to use the laws of large numbers in
order to replace $(A(k,t))_{k\geq 1}$ by more manipulable quantities and
propose to investigate the error made during this approximation. The first
possible approximation is the following.

Approximation 1:

$(A(k,t))_{k\geq 1}\approx(c_{k})_{\geq 1}e^{\alpha t}\mathcal{E}.$

However, this is unsatisfactory for practical applications since the random
variable $\mathcal{E}$ is not observable at finite times. Another idea is to
exploit the fact that the random variable appearing in the law of large
numbers for $A(k,t)$ and for $N_{t}$ is the same. This leads to the second
approximation.

Approximation 2:

$A(k,t)\approx(c_{k})_{\geq 1}N_{t}.$

In order to investigate the errors made during this approximation (at least
asymptotically), one would like to have central limit theorems associated to
the law of large numbers for the frequency spectrum. In a previous work [13],
we have showed that the error in the convergence of $e^{-\alpha t}N_{t}$ is of
order $e^{\alpha t/2}$ and obtained a central limit theorem for this error. An
important aspect of the method introduced in [13] is that it can be used to
derive CLTs for other branching processes counted by random characteristics.
In particular, the main goal of this work is to obtain central limit theorems
for the convergence of the frequency spectrum. We also study the Markovian
cases (when $\mathbb{P}_{V}$ is exponential) where we can obtain more
information on the limit distribution.

The original motivation of this study (and of other works on this model [6, 7,
5]) comes from the works of Sabeti and al. [18] where the frequency spectrum
is used to detecte positive selection of an allele in an increasing
population. More specifically, suppose that you want to detect the positive
selection of an allele on a given gene. The main idea is that, under neutral
evolution, the allele under consideration needs a long time to reach a high
frequency in the population. Hence, if the frequency of the allele w.r.t. its
age is significantly higher than the expected frequency (w.r.t. its age and
under neutral growth), this anomaly would suggest a positive selection of this
allele. The main problem is now to be able to estimate how old the allele is.
Sabeti and al. remarked that the allelic partition can be used as a clock to
estimate the age of an allele. More precisely, their study begins by selecting
a small region of chromosome which characterized the presence of the allele
under consideration. Now, the type of an individual, at a distance $x$
(measured in kb) from the core region, is the sequence of $x$ bases following
the core region (excluded). As a consequence, the allelic partition of the
subpopulation carrying the allele becomes thinner as $x$ increases (because
the higher $x$ is, the higher is the probability that a mutation occurred on
the sequence of $x$ bases). Finally, the speed of fragmentation of the allelic
partition, when $x$ increases, gives clues on the age of the allele. One of
the purposes of this model is to understand how the frequency spectrum evolves
under neutral evolution. In this work, we discuss some aspects of this method
and give some directions in order to construct rigorous tests for the positive
selection (see Section 7).

The paper is organized as follows. Section 2 is devoted to the mathematical
description of the model and to preliminary results which are used in the
sequel. Section 3 gives the mains theoretical results of this work and, in
particular, a central limit theorems which allow to study the error in our
proposed approximations. Section 4, 5, 6 are devoted to the proofs of Theorem
3.1, 3.3 and 3.5 respectively. Finally, in Section 7 we perform some numerical
studies on the model to stress the quality of our approximation. The
discussions about the method of Sabeti and al. are given in this last section.
An appendix contains some technical proofs and a section which is a reminder
of renewal theory.

##  2 Model and preliminaries

In this work, we consider a branching population with the following dynamic:
starting with a single individual (called the _ancestor_) whose lifetime is
distributed according to an arbitrary probability distribution
$\mathbb{P}_{V}$, this _ancestor_ gives birth to new individuals at a
Poissonian rate $b$. Each birth event giving a single new individual. From
this point, each child of the ancestor lives and gives birth according to the
same mechanism independently from the other individuals in the population.
This formal description can be made rigorous through the definition of a
probability distribution on the set of chronological trees. For the details of
such construction, we refer the reader to [14]. The first quantity of interest
when studying such population is the number $N_{t}$ of alive individuals in
the population at a fixed time $t$ (assuming that the time $t=0$ is birth-date
of the ancestor). The process $(N_{t},\ t\in\mathbb{R}_{+})$ is known as
binary homogeneous Crump-Mode-Jagers process and is a simple example of non-
Markovian branching process. In the sequel, we denote by $W(t)$ the
expectation of $N_{t}$ conditionally on the non-extinction at time $t$. That
is

$W(t):=\mathbb{E}\left[N_{t}\mid N_{t}>0\right].$

In [14], the author shows that the random variable $N_{t}$ is geometrically
distributed under $\mathbb{P}_{t}$ with parameter $\frac{1}{W(t)}$. In
addition, the author of [14] showed that the Laplace transform of $W$ can be
linked to the Laplace transform of $\mathbb{P}_{V}$ through the relation

$\int_{[0,\infty)}W(s)e^{-\lambda s}\ ds=\frac{1}{\psi(\lambda)},\
\forall\lambda>\alpha,$

where

$\psi(x)=x-\int_{(0,\infty]}\left(1-e^{-rx}\right)b\mathbb{P}_{V}(dr),\ \
x\in\mathbb{R}_{+},$ (2.1)

and $\alpha$ is the largest root of $\psi$. In particular, the Laplace
transform of $\mathbb{P}_{V}$ can be expressed in terms of $\psi$,

$\int_{\mathbb{R}_{+}}e^{-\lambda
v}\mathbb{P}_{V}(dv)=1+\frac{\psi(\lambda)-\lambda}{b}.$ (2.2)

In this work, we assume that $\alpha$ is a strictly positive real number. This
case is called the supercritical case and is equivalent to $b\mathbb{E}[V]>1$.
In the supercritical case, the real number $\alpha$ is called the Malthusian
parameter of the population because it corresponds to the mean exponential
growth rate of the population. Before gong further, let us remark that
equation (2.2) leads easily to the following identity:

$\int_{\mathbb{R}_{+}}e^{-\alpha v}\mathbb{P}_{V}(dv)=1-\frac{\alpha}{b}.$
(2.3)

Many previous works demonstrate [6, 7, 8] that some properties of the
splitting tree were easier to study on the tree describing only the
genealogical relation between the lineages of the individuals alive at time
$t$. For instance, in the model with mutations, the difference between two
individuals in term of type lies only on the time past since their lineages
has diverged. Hence, this particular genealogical tree, known as _coalescent
point processes_ (CPP), contains the essential information to study the
allelic partition. In order to derive the law of that genealogical tree, one
needs to characterize the joint law of the _times of coalescence_ between
pairs of individuals in the population, which are the times since their
lineages have split.

In [14], the author defines an order on the set of individuals alive at a
fixed time $t$ and consider the sequence of times of coalescences
$(H_{i})_{0\leq i\leq N_{t}-1}$ between two consecutive individuals (that is
$H_{i}$ is the time passed since the lineage of individuals $i$ and $i+1$ have
diverged) with the convention that the older lineage is the first one (i.e.
$H_{0}=t$). Moreover, in [14], the author shows that the random vector
$(H_{i})_{0\leq i\leq N_{t}-1}$ can be produced from a sequence
$(H_{i})_{i\geq 1}$ of i.i.d. random variable stopped at its first value
greater than $t$ and such that

$\mathbb{P}\left(H_{1}>s\right)=\frac{1}{W(s)},\quad s\in\mathbb{R}_{+}.$

To summarize, given the population is still alive at time $t$, one can forget
about the details of the splitting tree and code the genealogy by a new object
called the coalescent point process (CPP). Its law is the law of a sequence
$\left(H_{i}\right)_{0\leq i\leq N_{t}-1}$, where the family
$\left(H_{i}\right)_{i\geq 1}$ is i.i.d. with the same law as $H$, stopped
before its first value $H_{N_{t}}$ greater than $t$, and $H_{0}$ is
deterministic equal to $t$ (see Figure 1).

Although we do not use directly the CPP in this work, this object allowed us
to obtain [5] formulas for the moments of the frequency spectrum which are
widely used in the sequel.

######  Remark 2.1.

Let $N$ be a integer valued random variable. In the sequel we said that a
random vector with random size $\left(X_{i}\right)_{1\leq i\leq N}$ form an
i.i.d. family of random variables independent of $N$, if and only if

$\left(X_{1},\dots,X_{N}\right)\overset{d}{=}\left(\tilde{X_{1}},\dots,\tilde{X}_{N}\right),$

where $\left(\tilde{X}_{i}\right)_{i\geq 1}$ is a sequence of i.i.d. random
variables distributed as $X_{1}$ independent of $N$.

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$12$$13$$14$$15$ Figure 1:  A coalescent
point process for $16$ individuals, hence $15$ branches. (Image by A. Lambert)

Before going further, let us point out that if we define $N_{t}$ as the first
value of the sequence $\left(H_{i}\right)_{i\geq 1}$ greater than $t$, i.e.

$N_{t}=\inf\\{i\geq 1\mid H_{i}>t\\},$

then $N_{t}$ is indeed geometric with the expected parameter. More precisely,
for a positive integer $k$,

$\mathbb{P}\left(N_{t}=k\mid
N_{t}>0\right)=\frac{1}{W(t)}\left(1-\frac{1}{W(t)}\right)^{k-1}.$ (2.4)

In particular,

$\mathbb{E}\left[N_{t}\mid N_{t}>0\right]=W(t).$ (2.5)

Moreover, it can be showed (see [17]), that

$\mathbb{E}N_{t}=W(t)-W\star\mathbb{P}_{V}(t),$ (2.6)

and

$\mathbb{P}\left(N_{t}>0\right)=1-\frac{W\star\mathbb{P}_{V}(t)}{W(t)},$ (2.7)

where

$W\star\mathbb{P}_{V}(t):=\int_{[0,t]}W(t-s)\mathbb{P}_{V}(ds).$

Now, let us introduce the mathematical formalism for the mutation process used
in this work (this formalism comes from [5]). Since only the mutations
occurring on the lineages of living individuals at time $t$ can be observed,
it follows from standard properties on Poisson point processes, that the
mutation process can be defined directly on the CPP. So, let $\mathcal{P}$ be
a Poisson random measure on $[0,t]\times\mathbb{N}$ with intensity measure
$\theta\lambda\otimes C$, where $C$ is the counting measure on $\mathbb{N}$,
then the mutation random measure $\mathcal{N}$ on the CPP is defined by

$\mathcal{N}\left(da,di\right)=\mathds{1}_{H_{i}>t-a}\mathds{1}_{i<\mathcal{N}_{t}}\mathcal{P}\left(di,da\right),$

where an atom at $(a,i)$ means that the $i$th branch experiences a mutation at
time $t-a$. We suppose that each individual inherits the type of its parent.
This rule yields a partition of the population by types. The distribution of
the sizes of the families in the population is called the frequency spectrum
and is defined as the sequence $\left(A(k,t)\right)_{k\geq 1}$ where $A(k,t)$
is the number of types carried by exactly $k$ individuals in the alive
population at time $t$, excluding the family holding the ancestral type of the
population (i.e. individuals holding the same type as the root at time $0$).
This last family is called _clonal_ , as the ancestral type.

In the study of the frequency spectrum, an important role is played by the law
of the clonal family. We denote by $Z_{0}(t)$ the size of this family at time
$t$.

To study this family, it is easier to consider the clonal splitting tree
constructed from the original splitting tree by cutting every branches beyond
mutations. This clonal splitting tree is a standard splitting tree without
mutations, where individuals are killed as soon as they die or experience a
mutation. The new lifespan law is therefore the minimum between an exponential
random variable of parameter $\theta$ and an independent copy of $V$. It is
straightforward by simple manipulations of Laplace transforms that the Laplace
exponent of the corresponding contour process is

$\psi_{\theta}(x)=x-\int_{(0,\infty]}\left(1-e^{-rx}\right)\Lambda_{\theta}(dr)=\frac{x\psi(x+\theta)}{x+\theta}.$

We denote by $W_{\theta}$ the corresponding scale function. This leads to,

$\mathbb{P}\left(Z_{0}(t)=k\mid
Z_{0}(t)>0\right)=\frac{1}{W_{\theta}(t)}\left(1-\frac{1}{W_{\theta}(t)}\right)^{k-1}.$

When $\alpha>\theta$ (resp. $\alpha=\theta$, $\alpha<\theta$), this new tree
is supercritical (resp. critical, sub-critical) and we talk about _clonal
supercritical case_ (resp. _critical_ , _sub-critical_ case).

Moreover, the law of $Z_{0}$ conditionally on the event $\\{N_{t}>0\\}$ can be
obtained, and is given by

$\mathbb{P}\left(Z_{0}(t)=k\mid N_{t}>0\right)=\frac{e^{-\theta
t}W(t)}{W_{\theta}(t)^{2}}\left(1-\frac{1}{W_{\theta}(t)}\right)^{k-1},\quad\forall
k\geq 1.$ (2.8)

For the rest of this paper, unless otherwise stated, the notation
$\mathbb{P}_{t}$ refers to $\mathbb{P}\left(.\mid N_{t}>0\right)$ whereas
$\mathbb{P}_{\infty}$ refers to the probability measure conditioned on the
non-extinction event (which has positive probability in the supercritical
case).

Finally, we recall the asymptotic behavior of the scale functions $W(t)$ and
$W_{\theta}(t)$ which is widely used in the sequel,

######  Lemma 2.2.

([6, Thm. 3.21]) There exist a positive constant $\gamma$ such that,

$e^{-\alpha t}\psi^{\prime}(\alpha)W(t)-1=\mathcal{O}\left(e^{-\gamma
t}\right).$

In the case that $\theta<\alpha$ (clonal supercritical case),

$W_{\theta}(t)\underset{t\to\infty}{\sim}\frac{e^{\left(\alpha-\theta\right)t}}{\psi_{\theta}(\alpha-\theta)}.$

In the case that $\theta>\alpha$ (clonal sub-critical case),

$W_{\theta}(t)=\frac{\theta}{\psi(\theta)}+\mathcal{O}\left(e^{-\left(\theta-\alpha\right)t}\right).$

In the case where $\theta=\alpha$ (clonal critical case),

$W_{\theta}(t)\underset{t\to\infty}{\sim}\frac{\theta
t}{\psi^{\prime}(\alpha)}.$

For a purpose, a more precise description of the asymptotic behavior of $W$ is
needed. It is given by the following result.

######  Lemma 2.3.

[13, Prop. 5.1] There exists a positive non-increasing càdlàg function $F$
such that

$W(t)=\frac{e^{\alpha t}}{\psi^{\prime}(\alpha)}-e^{\alpha t}F(t),\quad t\geq
0,$

and

$\lim\limits_{t\to\infty}e^{\alpha
t}F(t)=\begin{cases}\frac{1}{b\mathbb{E}V-1}&\mbox{if}\ \mathbb{E}V<\infty,\\\
0&\mbox{otherwise.}\end{cases}$

From this Lemma and (2.7), one can easily deduce that

$\mathbb{P}\left(\text{NonEx}\right)=\lim\limits_{t\to\infty}\mathbb{P}\left(N_{t}>0\right)=\frac{\alpha}{b},$
(2.9)

where NonEx refer to the non-extinction event.

In [5], we show that a CPP stopped at time $t$ with scale function $W$ can be
constructed by grafting independent CPP stopped at a fixed time $a\leq t$ on a
CPP stopped at time $t-a$ with an explicit scale function different of $W$
(see Figure 2).

tt-a$\mathcal{P}^{(1)}$$\mathcal{P}^{(2)}$$\mathcal{P}^{(3)}$$\mathcal{P}^{(4)}$
Figure 2:  Adjunction of trees.

Moreover, we showed that the frequency spectrum can be expressed as an
integral with respect to the random measure $\mathcal{N}$ along the CPP, that
is

$\prod_{i=1}^{l}A(k_{i},t)=\sum_{i=1}^{l}\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{(u)}_{0}(a)=k_{i}}\sum_{u_{1:l-1}=1}^{N^{(t)}_{t-a}}\prod_{\underset{i\neq
j}{j=1}}^{l-1}A^{(u_{j})}(k_{j},a)\ \mathcal{N}\left(da,du\right),$ (2.10)

where $A^{(u)}(k,a)$ (resp. $Z^{(u)}_{0}$) refers to the frequency spectrum
(resp. clonal family) of the $u$th grafted sub-CPP, and
$\sum_{u_{1:l-1}=1}^{N^{(t)}_{t-a}}$ denotes for the multi-sum

$\sum_{u_{1}=1}^{N^{(t)}_{t-a}}\dots\sum_{u_{l-1}=1}^{N^{(t)}_{t-a}}.$

Moreover, in [5, Thm, 3.1] we show that the expectation of such integral can
be computed easily when the integrand presents local independence properties
with the random measure as in formula (2.10). Equation (2.10) is used later to
obtain some moments estimates useful to prove our theorems. In particular,
this allows to prove that (see [5]) for any positive integer $k$ and $l$,

$\mathbb{E}_{t}A(k,t)=W(t)\int_{0}^{t}\frac{\theta e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}ds,$ (2.11)

and

$\displaystyle\mathbb{E}A(k,t)A(l,t)$
$\displaystyle=2W(t)^{2}\int_{0}^{t}\frac{\theta e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}ds\int_{0}^{t}\frac{\theta
e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{l-1}ds$
$\displaystyle-W(t)\int_{0}^{t}2\theta\ \frac{e^{-\theta
a}W(a)}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{l-1}\int_{0}^{s}\frac{\theta
e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}dsda$
$\displaystyle-W(t)\int_{0}^{t}2\theta\ \frac{e^{-\theta
a}W(a)}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\int_{0}^{s}\frac{\theta
e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{l-1}dsda$
$\displaystyle+W(t)\mathbb{E}\int_{0}^{t}\theta
W(a)^{-1}\left(\mathbb{E}\left[A(k,t)\mathds{1}_{Z_{0}(a)=l}\right]+\mathbb{E}\left[A(l,t)\mathds{1}_{Z_{0}(a)=k}\right]\right)da$
$\displaystyle+\mathds{1}_{l=k}W(t)\int_{0}^{t}\frac{\theta e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}ds.$ (2.12)

These tools also allow, for instance, to prove next two results [14, 8, 5].

######  Theorem 2.4.

There exists a random variable $\mathcal{E}$, such that

$\lim\limits_{t\to\infty}e^{-\alpha
t}N_{t}=\frac{\mathcal{E}}{\psi^{\prime}(\alpha)},\quad a.s.\text{ and in
}L^{2}.$

Moreover, under $\mathbb{P}_{\infty}$, $\mathcal{E}$ is exponentially
distributed with parameter one.

######  Theorem 2.5.

For any positive integer $k$,

$\lim\limits_{t\to\infty}e^{-\alpha
t}A(k,t)=\frac{c_{k}\mathcal{E}}{\psi^{\prime}(\alpha)},\quad a.s.\text{ and
in }L^{2},$

where $\mathcal{E}$ is the random variable of the Theorem 2.4 and

$c_{k}=\int_{0}^{\infty}\frac{\theta e^{-\theta
a}}{W_{\theta}(a)}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}da.$ (2.13)

##  3 Main results

The a.s. convergence stated in Section 2 suggests studying the second order
properties of the convergence to get central limit theorem. Our main result,
Theorem 3.5, allows to study the asymptotic error in the approximation $2$
proposed in the introduction of this work. In addition, we prove more standard
central limit theorems which are interesting from the theoretical point of
view.

Before going further, we recall that the Laplace distribution with mean
$\mu\in\mathbb{R}^{n}$ and covariance matrix $K$ is the probability
distribution whose characteristic function is given, for all
$\lambda\in\mathbb{R}^{n}$ by

$\frac{1}{1+\frac{1}{2}\lambda^{\prime}K\lambda-i\mu^{\prime}\lambda}$

We denote this law by $\mathcal{L}\left(\mu,K\right)$. We also recall that, if
$G$ is a Gaussian random vector with mean $\mu$ and covariance matrix $K$ and
$\mathcal{E}$ is an exponential random variable with parameter $1$ independent
of $G$, then $\sqrt{\mathcal{E}}G$ is Laplace $\mathcal{L}\left(\mu,K\right)$.

###  3.1 CLT for the convergence of Theorem 2.5

######  Theorem 3.1.

Suppose that $\theta>\alpha$ and
$\int_{[0,\infty)}e^{\left(\theta-\alpha\right)v}\mathbb{P}_{V}(dv)>1$ . Then,
we have, under $\mathbb{P}_{\infty}$,

$\left(e^{\alpha\frac{t}{2}}\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\right)_{k\in\mathbb{N}}\xrightarrow[t\to\infty]{(d)}\mathcal{L}\left(0,K\right),$

where $K$ is some covariance matrix and the constants $c_{k}$ are defined in
(2.13).

The proof of this result can be found in Section 4.

######  Remark 3.2.

We are not able to compute explicitly the covariance matrix $K$ in the general
case due to our method of demonstration. However, all our other results give
explicit formulas. In particular, the case where $\mathbb{P}_{V}$ is
exponential is given by the next theorem. The Yule case is also covered in the
following theorem for $d=0$ although it does not satisfy the hypothesis of
Theorem 3.1.

######  Theorem 3.3.

Suppose that $V$ is exponentially distributed with parameter $d\in[0,b)$. In
this case, $\alpha=b-d$. We still suppose that $\alpha<\theta$, then

$\left(e^{\alpha\frac{t}{2}}\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\right)_{k\in
N}\xrightarrow[t\to\infty]{(d)}\mathcal{L}\left(0,K\right),\ \text{w.r.t.}\
\mathbb{P}_{\infty},$

where $K$ is given by

$K_{l,k}=M_{l,k}+c_{k}c_{l}\frac{\alpha}{b}\left(1-6\frac{d}{\alpha}\right),$

and

$\displaystyle M_{l,k}=$ $\displaystyle
2\psi^{\prime}(\alpha)\int_{0}^{\infty}\frac{\theta e^{-\theta
a}}{W_{\theta}(a)^{2}}\left(\left(1-\frac{1}{W_{\theta}(a)}\right)^{l-1}\left(\mathbb{E}_{a}\left[A(k,a)\right]-c_{k}W(a)\right)+\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\left(\mathbb{E}_{a}\left[A(l,a)\right]-c_{l}W(a)\right)\right)da$
$\displaystyle-\psi^{\prime}(\alpha)\int_{0}^{\infty}\theta
W(a)^{-1}\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\mathds{1}_{Z_{0}(a)=l}+\left(A(l,a)-c_{l}N_{a}\right)\mathds{1}_{Z_{0}(a)=k}\right]$
$\displaystyle+\mathds{1}_{l=k}\int_{0}^{\infty}\frac{\theta e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}ds,$ (3.1)

where $W$, $W_{\theta}$, $\psi^{\prime}(\alpha)$ are defined in the Section 2.

The proof of this result can be found in Section 6. Note that an explicit
formula for $\mathbb{E}_{t}A(k,t)$ is given by (2.11). Explicit formulas for
$\mathbb{E}_{t}\left[A(k,t)\mathds{1}_{Z_{0}(t)=l}\right]$ can also be found
in Proposition 4.5 of [5], and a formula for
$\mathbb{E}_{t}\left[N_{a}\mathds{1}_{Z_{0}(t)=k}\right]$ can be found in
Proposition 4.1 of [6].

######  Remark 3.4.

The condition on $V$ in Theorem 3.1 is required only to ensure controls of the
moments of the considered quantities. However, although the Yule case does not
satisfy this condition ($V=\infty$ p.s.) it is included in this last theorem
(d=0). This suggests that the condition on $V$ may not be needed.

###  3.2 CLT for the error between $A(k,t)$ and $c_{k}N_{t}$

The next theorem concerns the error between $A(k,t)$ and $c_{k}N_{t}$. Once
again, we have an explicit expression of the covariance matrix of the limit.

######  Theorem 3.5.

Suppose that $\theta>\alpha$, then

$\psi^{\prime}(\alpha)\left(e^{\alpha\frac{t}{2}}\left(A(k,t)-c_{k}N_{t}\right)\right)_{k\in\mathbb{N}}\xrightarrow[t\to\infty]{(d)}\mathcal{L}\left(0,M\right),\
\text{w.r.t.}\ \mathbb{P}_{\infty},$

where $M$ is defined in relation (3.3).

The proof of this result can be found in Section 5.

######  Remark 3.6.

We do not known yet if the exponential random variable appearing the Gaussian
mixing leading to a Laplace distribution is the same as the exponential limit
of $e^{-\alpha t}A(k,t)$. However, the CLT for Markov branching processes in
[3] suggest that it is, actually, the case. If, this is true in our case, it
would be enough to know the correlations between the limits involved in
Theorem 3.1 and 3.5 to obtain an explicit expression for the covariance matrix
in Theorem 3.1.

##  4 Proof of Theorem 3.1

The proof of this theorem is based on the proof of the central limit theorem
for the process $(N_{t},\ t\in\mathbb{R}_{+})$ given in [13]. The structure of
the proof follows the same lines and is detailed in Section 4 of [13]. In a
sake of conciseness, we only highlight the difficulties arising in our new
context. The results which are straightforward rewording of the proofs given
in [13] are left to the reader. However, we think it is necessary to recall
some aspects of [13], in particular from [13, Section 4]. First, we recall
that there exists a family $(N^{(i)}_{t},\ t\in\mathbb{R}_{+})_{i\geq 1}$ of
i.i.d. population counting processes with the same law as $\left(N_{t},\
t\in\mathbb{R}_{+}\right)$, and a Poisson random measure $\xi$ on
$\mathbb{R}_{+}$ with intensity $b\,da$ such that

$N_{t}=\int_{[0,t]}N^{(\xi_{u})}_{t-u}\mathds{1}_{V_{\emptyset}>u}\
\xi(du)+\mathds{1}_{V_{\emptyset}>t},\quad\text{almost surely},$ (4.1)

where $\xi_{u}=\xi\left([0,u]\right)$. In addition, we have that
$t\to\mathbb{E}\left[N_{t}\mathcal{E}\right]$ is the unique solution bounded
on finite intervals of the renewal equation,

$\displaystyle f(t)=$ $\displaystyle\int_{\mathbb{R}_{+}}f(t-u)be^{-\alpha
u}\mathbb{P}\left(V>u\right)du$ $\displaystyle+\alpha
b\mathbb{E}\left[N_{\cdot}\right]\star\left(\int_{\mathbb{R}_{+}}e^{-\alpha
v}\mathbb{P}\left(V>\cdot,V>v\right)dv\right)(t)$
$\displaystyle+\alpha\int_{\mathbb{R}_{+}}e^{-\alpha
v}\mathbb{P}\left(V>t,V>v\right)dv,$ (4.2)

and it is given by

$\mathbb{E}\left[N_{t}\mathcal{E}\right]=\left(1+\frac{\alpha}{b}-e^{-\alpha
t}\right)W(t)-\left(1-e^{-\alpha t}\right)W\star\mathbb{P}_{V}(t).$ (4.3)

We also recall that equation (4.2) is obtained by taking the product
$N_{t}N_{s}$, for some real number $t$ and $s$. Now, equation (4.1) allows to
obtain a renewal equation for $\mathbb{E}[N_{t}N_{s}]$ which leads to (4.2)
when taking the limit in $s$ of the renormalized equation. We also recall that
Lemma 2.3 and equation (2.7) gives

$\frac{1}{\mathbb{P}\left(N_{t}>0\right)}=\frac{b}{\alpha}-\frac{b\mu\psi^{\prime}(\alpha)}{\alpha}e^{-\alpha
t}+o(e^{-\alpha t}).$ (4.4)

This also leads, in conjunction with equation (4.3), to

$\mathbb{E}_{t}N_{t}\mathcal{E}=\frac{2e^{\alpha
t}}{\psi^{\prime}(\alpha)}-\frac{1}{\psi^{\prime}(\alpha)}-3\mu+o(1).$ (4.5)

Finally, let us recall that for any fixed time $u$, there is a natural order
(for instance given by the contour process [14]) of the individuals alive at
this time. Moreover, we denote, for $1\leq i\leq N_{t}$, $O_{i}^{(u)}$ the
residual lifetime of the $i$th individual alive at time $u$. The law of the
vector $(O^{(u)}_{2},\dots,O^{(u)}_{N_{u}})$ is given by the following lemma
which comes from [13].

######  Lemma 4.1.

Let $u$ in $\mathbb{R}_{+}$, we denote by $O_{i}$ for $i$ an integer between
$1$ and $N_{u}$ the residual lifetime of the $i$th individuals alive at time
$u$. Then under $\mathbb{P}_{u}$, the family $\left(O_{i},\
i\in\\{1,\dots,N_{u}\\}\right)$ form a family of independent random variables,
independent of $N_{u}$, and, expect $O_{1}$, having the same distribution,
given by, for $2\leq i\leq N_{t}$,

$\mathbb{P}_{u}(O_{i}\in dx)=\int_{\mathbb{R}_{+}}\
\frac{W(u-y)}{W(u)-1}b\mathbb{P}\left(V-y\in dx\right)\ dy.$ (4.6)

Moreover, it follows that the family
$\left(N_{s}(O_{i}),s\in\mathbb{R}_{+}\right)_{1\leq i\leq N_{u}}$ is an
independent family of process, i.i.d. for $i\geq 2$, and independent of
$N_{u}$.

To end this reminder, let us recall the decomposition of the limiting random
variable $\mathcal{E}$ (given for instance in Theorem 2.4) at a fixed time
$u$.

######  Lemma 4.2.

[13, Lemma 6.8] We have the following decomposition of $\mathcal{E}$,

$\mathcal{E}=e^{-\alpha
u}\sum_{i=1}^{N_{u}}\mathcal{E}_{i}\left(O_{i}\right),\quad a.s.$

Moreover, under $\mathbb{P}_{u}$, the random variables
$\left(\mathcal{E}_{i}\left(O_{i}\right)\right)_{i\geq 1}$ are independent,
independent of $N_{u}$, and identically distributed for $i\geq 2$.

We can now start the proof of theorem 3.1. As in [13], the proof begins by
some estimate on moments.

###  4.1 Preliminary moments estimates

We start by computing the moment in the case of a standard splitting tree.
According to [13, Section 4], the next step is to obtain the same kind of
estimates in the case of a splitting tree whose ancestor individual has a
lifetime distribution which can be different from the rest of the population.

####  4.1.1 Case $V_{\emptyset}\overset{\mathcal{L}}{=}V$

One of the main difficulties to extend the preceding proof to the frequency
spectrum is to get estimates on

$\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)^{n}\right],\ \text{for }n=2\text{ or }3.$

We first study the renewal equation satisfied by $\mathbb{E}A(k,t)\mathcal{E}$
similarly as in [13, Lemma 6.1].

######  Lemma 4.3 (Joint moment of $\mathcal{E}$ and $A(k,t)$).

$\mathbb{E}\left[A(k,t)\mathcal{E}\right]$ is the unique solution bounded on
finite intervals of the renewal equation,

$\displaystyle f(t)=$ $\displaystyle\int_{\mathbb{R}_{+}}f(t-u)be^{-\alpha
u}\mathbb{P}\left(V>u\right)du$
$\displaystyle+\alpha\mathbb{E}\left[A(k,.)\right]\star
b\left(\int_{\mathbb{R}_{+}}e^{-\alpha
v}\mathbb{P}\left(V>.,V>v\right)dv\right)(t)$
$\displaystyle+\alpha\mathbb{E}\left[\mathcal{E}X_{t}\right],$ (4.7)

with $X_{t}$ the number of families of size $k$ alive at time $t$ whose
original mutation has taken place during the lifetime of the ancestor
individual.

###### Proof.

We recall that $A(k,t)$ is the number of non-ancestral families of size $k$ at
time $t$. Similarly, as for $N_{t}$, $A(k,t)$ can be obtained as the sum of
the contributions of all the trees grafted on the lifetime of the ancestor
individual in addition to the mutations which take place on the ancestral
branch, that is,

$A(k,t)=\int_{[0,t]}A(k,t-u,\xi_{u})\mathds{1}_{V_{\emptyset}>u}\xi(du)+X_{t},$

where $\left(A(k,t,i),t\in\mathbb{R}_{+}\right)_{i\geq 1}$ is a family of
independent processes having the same law as $A(k,t)$. Now, taking the product
$A(k,t)N_{s}$ and using the same arguments as in the proof of lemma [13, Lemma
6.1] to take the limit in $s$ leads to the result. In particular, the last
term is obtained using that

$\lim\limits_{s\to\infty}\mathbb{E}\left[X_{t}\frac{N_{s}}{W(s)}\right]=\mathbb{E}\left[X_{t}\mathcal{E}\right].$

∎

The result of Lemma 4.3 is quite disappointing since the presence of the
mysterious process $X_{t}$ prevents any explicit resolution of equation (4.3).
However, one may note that equation (4.3) is quite similar to equation (4.2)
driving $\mathbb{E}N_{t}\mathcal{E}$, so if the contribution of $X_{t}$ in the
renewal structure of the process is small enough, one can expect the same
asymptotic behavior for $\mathbb{E}A(k,t)\mathcal{E}$ as for
$\mathbb{E}N_{t}\mathcal{E}$. Moreover, we clearly have on $X_{t}$ the
following a.s. estimate,

$X_{t}\leq\int_{[0,t]}\mathds{1}_{Z^{(u)}_{0}(t-u)>0}\mathds{1}_{V>u}\xi(du),$
(4.8)

where $Z^{(i)}_{0}$ denote for the ancestral families on the $i$th trees
grafted on the ancestral branch. Hence, if we take $\theta>\alpha$ and we
suppose $V<\infty$ a.s., one can expect that $X_{t}$ decreases very fast.
These are the ideas the following Lemma is based on. Moreover, as it is seen
in the proof of the following lemma, the hypothesis $V<\infty$ a.s. can be
weakened.

######  Lemma 4.4.

Under the hypothesis of Theorem 3.1, for all $k\geq 1$, there exists a
constant $\gamma_{k}\in\mathbb{R}$ such that,

$\lim\limits_{t\to\infty}\mathbb{E}N_{t}\mathcal{E}c_{k}-\mathbb{E}A(k,t)\mathcal{E}=\gamma_{k}.$
(4.9)

###### Proof.

Combining equations $\eqref{eq:renewNt}$ and $\eqref{eq:renewAk}$, we get
that,

$\displaystyle\mathbb{E}N_{t}\mathcal{E}c_{k}-\mathbb{E}A(k,t)\mathcal{E}=\int_{\mathbb{R}_{+}}\left(\mathbb{E}N_{t-u}\mathcal{E}c_{k}-\mathbb{E}A(k,t-u)\mathcal{E}\right)be^{-\alpha
u}\mathbb{P}\left(V>u\right)du$ $\displaystyle\quad+\underbrace{\alpha
b\left(c_{k}\mathbb{E}N_{.}-\mathbb{E}\left[A(k,.)\right]\right)\star\left(\int_{\mathbb{R}_{+}}e^{-\alpha
v}\mathbb{P}\left(V>.,V>v\right)dv\right)(t)}_{:=\xi^{(k)}_{1}(t)}$
$\displaystyle\quad\quad+\underbrace{c_{k}\mathbb{P}\left(V>t\right)-\alpha\mathbb{E}\left[X_{t}\mathcal{E}\right]}_{:=\xi^{(k)}_{2}(t)},$

which is also a renewal equation. On one hand, using equations (2.5) and
(2.11) imply that

$\mathbb{E}_{t}\left[c_{k}N_{t}-A(k,t)\right]=W(t)\int_{t}^{\infty}\frac{\theta
e^{-\theta s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}\
ds,$

which leads using Lemma 2.2, to

$\displaystyle\xi_{1}(t)=$
$\displaystyle\alpha\int_{\mathbb{R}_{+}}\left(c_{k}\mathbb{E}N_{t-u}-\mathbb{E}\left[A(k,t-u)\right]\right)\int_{\mathbb{R}_{+}}e^{-\alpha
v}\mathbb{P}\left(V>u,V>v\right)dvdu$
$\displaystyle\leq\mathcal{C}\int_{[0,t]}e^{(\alpha-\theta)t-u}\mathbb{P}\left(V>u\right)du\int_{[0,\infty)}e^{-\alpha
u}du$
$\displaystyle\leq\frac{\mathcal{C}}{\alpha}e^{-\left(\theta-\alpha\right)t}\int_{0}^{t}e^{(\theta-\alpha)u}\mathbb{P}\left(V>u\right)du,$
(4.10)

for some positive real constant $\mathcal{C}$.

The derivative of the r.h.s. of (4.1.1) is given by

$\frac{\mathcal{C}}{\alpha}e^{-(\theta-\alpha)t}\left(e^{(\theta-\alpha)t}\mathbb{P}\left(V>t\right)-(\alpha-\theta)\int_{0}^{t}e^{(\theta-\alpha)u}\mathbb{P}\left(V>u\right)du\right),\quad
t>0,$ (4.11)

which is equal to

$\frac{\mathcal{C}}{\alpha}e^{-(\theta-\alpha)t}\left(1-\int_{[0,t]}e^{(\theta-\alpha)s}\mathbb{P}_{V}(ds)\right),\quad
t>0,$

using Stieljes integration by parts. Now, since,

$\int_{[0,\infty)}e^{\left(\theta-\alpha\right)s}\mathbb{P}_{V}(ds)>1,$

this shows that the right hand side of (4.1.1) is decreasing for $t$ large
enough. Moreover, it is straightforward to shows that the r.h.s. of (4.1.1) is
also integrable. This implies that $\xi_{1}^{(k)}$ is DRI from the same Lemma.
On the other hand, it follows from (4.8) that

$X_{t}\mathcal{E}\leq\mathcal{E}\int_{[0,t]}\mathds{1}_{Z^{(u)}_{0}(t-u)>0}\mathds{1}_{V>t}\xi(du).$
(4.12)

Then, we obtain using Cauchy-Schwarz inequality, that

$\mathbb{E}\left[X_{t}\mathcal{E}\right]\leq\sqrt{\frac{2\alpha}{b}}\mathbb{E}\left[\left(\int_{[0,t]}\mathds{1}_{Z^{(u)}_{0}(t-u)>0}\mathds{1}_{V>t}\xi(du)\right)^{2}\right]^{1/2}.$

It follows that we need to investigate the behavior of

$\mathbb{E}\left[\left(\int_{(0,t)}\mathds{1}_{Z^{(u)}_{0}(t-u)>0}\mathds{1}_{V>t}\xi(du)\right)^{2}\right],$

which is equal to

$\int_{0}^{t}\mathbb{P}\left(Z_{0}(t-u)>0\right)\mathbb{P}\left(V>t\right)bdu+\int_{[0,t]^{2}}\mathbb{P}\left(Z_{0}(t-v)>0\right)\mathbb{P}\left(Z_{0}(t-u)>0\right)\mathbb{P}\left(V>u,V>v\right)b^{2}du\
dv,$

using [13, Lemma 2.6]. Then, since, from (2.8) and Lemma 2.2,

$\mathbb{P}_{t-u}\left(Z_{0}(t-u)>0\right)=\frac{e^{-\theta(t-u)}W(t-u)}{W_{\theta}(t-u)}=\mathcal{O}(e^{-(\theta-\alpha)(t-u)}),$

it follows, using that the right hand side of (4.1.1) is DRI and Lemma A.1,
that $\xi^{(k)}_{2}$ is DRI. Finally, it comes from Theorem A.2, that

$\lim\limits_{t\to\infty}\mathbb{E}N_{t}\mathcal{E}c_{k}-\mathbb{E}A(k,t)\mathcal{E}=\frac{\alpha}{\psi^{\prime}(\alpha)}\int_{\mathbb{R}_{+}}\xi_{1}^{(k)}(s)+\xi^{(k)}_{2}(s)ds.$
(4.13)

∎

Using the preceding lemma, we can now get the quadratic error in the
convergence of the frequency spectrum.

######  Lemma 4.5 (Quadratic error for the convergence of $A(k,t)$.).

Let $k$ and $l$ two positive integers. Then under the hypothesis of Theorem
3.1, there exists a family of real numbers $\left(a_{k,l}\right)_{l,k\geq 1}$
such that,

$\lim\limits_{t\to\infty}e^{-\alpha
t}\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}\mathcal{E}c_{k}\right)\left(\psi^{\prime}(\alpha)A(l,t)-e^{\alpha
t}\mathcal{E}c_{l}\right)\right]=\frac{\alpha}{b}a_{k,l},$

where the sequence $\left(c_{k}\right)_{k\geq 1}$ is defined by (2.13).

###### Proof.

Now, noting

$c_{k}(t):=\int_{0}^{t}\frac{\theta e^{-\theta
a}}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}da,$ (4.14)

we have, from (2) and Lemma 2.3,

$\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}\left[A(k,t)A(l,t)\right]=2e^{2\alpha
t}c_{k}(t)c_{l}(t)+e^{\alpha t}\left(4\psi^{\prime}(\alpha)e^{\alpha
t}F(t)c_{k}(t)c_{l}(t)+\frac{R}{\psi^{\prime}(\alpha)}\right)+\mathcal{O}\left(1\right),$
(4.15)

with

$\displaystyle R:=$
$\displaystyle-\psi^{\prime}(\alpha)\int_{0}^{\infty}2\theta\ \frac{e^{-\theta
a}W(a)}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{l-1}\int_{0}^{a}\frac{e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}dsda$
$\displaystyle-\psi^{\prime}(\alpha)\int_{0}^{\infty}2\theta\ \frac{e^{-\theta
a}W(a)}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\int_{0}^{a}\frac{e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{l-1}dsda$
$\displaystyle+\psi^{\prime}(\alpha)\int_{0}^{\infty}\theta
W(a)^{-1}\left(\mathbb{E}_{t}\left[A(k,t)\mathds{1}_{Z_{0}(a)=l}\right]+\mathbb{E}_{t}\left[A(l,t)\mathds{1}_{Z_{0}(a)=k}\right]\right)da,$

and $F$, $\mu$ are defined in Lemma 2.3. Now, using (4.4), we have

$\mathbb{E}_{t}\mathcal{E}^{2}-2=-2\mu\psi^{\prime}(\alpha)e^{-\alpha
t}+o(e^{-\alpha t}),$

which leads to

$\displaystyle\mathbb{E}_{t}$ $\displaystyle\left[\left(e^{-\alpha
t}\psi^{\prime}(\alpha)A(k,t)-\mathcal{E}c_{k}\right)\left(e^{-\alpha
t}\psi^{\prime}(\alpha)A(l,t)-\mathcal{E}c_{l}\right)\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{t}\left[e^{-2\alpha
t}\psi^{\prime}(\alpha)^{2}A(k,t)A(l,t)\right]-c_{l}\mathbb{E}_{t}\left[e^{-\alpha
t}\psi^{\prime}(\alpha)A(k,t)\mathcal{E}\right]-c_{k}\mathbb{E}_{t}\left[e^{-\alpha
t}\psi^{\prime}(\alpha)A(l,t)\mathcal{E}\right]$
$\displaystyle+2c_{k}c_{l}-2c_{k}c_{l}\mu\psi^{\prime}(\alpha)e^{-\alpha
t}+o(e^{-\alpha t}),$ $\displaystyle=$ $\displaystyle
2\left(c_{k}(t)-c_{k}\right)\left(c_{l}(t)-c_{l}\right)-4\mu\psi^{\prime}(\alpha)c_{k}c_{l}e^{-\alpha
t}+Re^{-\alpha t}$
$\displaystyle-\left(2c_{k}(t)c_{l}+2c_{l}(t)c_{k}-2c_{k}c_{l}\psi^{\prime}(\alpha)e^{-\alpha
t}\mathbb{E}_{t}N_{t}\mathcal{E}\right)$
$\displaystyle+\psi^{\prime}(\alpha)c_{l}e^{-\alpha
t}\mathbb{E}_{t}\left[\left(c_{k}N_{t}-A(k,t)\right)\mathcal{E}\right]+\psi^{\prime}(\alpha)c_{k}e^{-\alpha
t}\mathbb{E}_{t}\left[\left(c_{l}N_{t}-A(l,t)\right)\mathcal{E}\right]+o(e^{-\alpha
t}),$

Since, by Lemma 2.2

$c_{k}(t)=c_{k}+\mathcal{O}(e^{-\theta t})=c_{k}+o(e^{-\alpha t}),$

it follows, combining (4.13), (4.5), and Lemma 4.4, that

$\displaystyle e^{\alpha t}\mathbb{E}_{t}$
$\displaystyle\left[\left(e^{-\alpha
t}\psi^{\prime}(\alpha)A(k,t)-\mathcal{E}c_{k}\right)\left(e^{-\alpha
t}\psi^{\prime}(\alpha)A(l,t)-\mathcal{E}c_{l}\right)\right]$ $\displaystyle=$
$\displaystyle\psi^{\prime}(\alpha)\left(c_{k}\gamma_{l}+c_{l}\gamma_{k}\right)+c_{k}c_{l}\left(2e^{\alpha
t}-2\psi^{\prime}(\alpha)\mathbb{E}_{t}N_{t}\mathcal{E}\right)+R-4\mu\psi^{\prime}(\alpha)c_{k}c_{l}+o(1)$
$\displaystyle=$
$\displaystyle\psi^{\prime}(\alpha)\left(c_{k}\gamma_{l}+c_{l}\gamma_{k}\right)+c_{k}c_{l}\left(\frac{1}{\psi^{\prime}(\alpha)}+3\mu\right)+R-4\mu\psi^{\prime}(\alpha)c_{k}c_{l}+o(1).$

The result follows readily from the fact that
$\mathbb{P}\left(N_{t}>0\right)\sim\frac{\alpha}{b}$.

∎

######  Lemma 4.6 (Boundedness of the third moment).

Let $k_{1},k_{2},k_{3}$ three positive integers, then

$\mathbb{E}\left[\prod_{i=1}^{3}\left|e^{-\frac{\alpha}{2}t}\left(\psi^{\prime}(\alpha)A(k_{i},t)-e^{\alpha
t}\mathcal{E}c_{k_{i}}\right)\right|\right]=\mathcal{O}\left(1\right).$

###### Proof.

We have,

$\mathbb{E}\left[\left|\prod_{i=1}^{3}\frac{\left(\psi^{\prime}(\alpha)A(k_{i},t)-e^{\alpha
t}\mathcal{E}c_{k_{i}}\right)}{e^{\frac{\alpha}{2}t}}\right|\right]\leq\prod_{i=1}^{3}\left(\mathbb{E}\left[\left|\frac{\left(\psi^{\prime}(\alpha)A(k_{i},t)-e^{\alpha
t}\mathcal{E}c_{k_{i}}\right)}{e^{\frac{\alpha}{2}t}}\right|^{3}\right]\right)^{\frac{1}{3}}.$

Hence, we only have to prove the Lemma for $k_{1}=k_{2}=k_{3}=k$. Hence,

$\displaystyle\mathbb{E}\left[\left|\frac{\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}\mathcal{E}c_{k}\right)}{e^{\frac{\alpha}{2}t}}\right|^{3}\right]$
$\displaystyle\leq
8\mathbb{E}\left[\left|\frac{\psi^{\prime}(\alpha)A(k,t)-c_{k}N_{t}}{e^{\frac{\alpha}{2}t}}\right|^{3}\right]+8c_{k}\mathbb{E}\left[\left|\frac{\psi^{\prime}(\alpha)N_{t}-N^{\infty}_{t}}{e^{\frac{\alpha}{2}t}}\right|^{3}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+8c_{k}\mathbb{E}\left[\left|\frac{N^{\infty}_{t}-e^{\alpha
t}\mathcal{E}}{e^{\frac{\alpha}{2}t}}\right|^{3}\right].$

The last two terms have been treated in the proof of [13, Lemma 6.4], and the
boundedness of

$\mathbb{E}\left[\left|\frac{\psi^{\prime}(\alpha)A(k,t)-c_{k}N_{t}}{e^{\frac{\alpha}{2}t}}\right|^{3}\right],$

follows from the following Lemma 4.7 and Hölder’s inequality.

∎

######  Lemma 4.7.

For all $k\geq 1$,

$\mathbb{E}\left[\left(\frac{A(k,t)-c_{k}N_{t}}{e^{-\frac{\alpha}{2}t}}\right)^{4}\right],$

is bounded.

Due to technicality, the proof of this lemma is postponed to the end in
appendix.

####  4.1.2 Arbitrary initial distribution case

The following Lemmas are the counter part of Lemmas 6.5, 6.6, and 6.7 of [13].
They play the same role in the proof of Theorem 3.1 as in the proof of the
central limit theorem given in [13]. In the sequel, we denote by
$\left(A(k,t,\Xi)\right)_{k\geq 1}$, the frequency spectrum of the splitting
tree where the lifetime of the ancestral individual is $\Xi$, in the same
manner as for $N_{t}\left(\Xi\right)$ in [13].

$\mathcal{E}_{i}:=\lim\limits_{t\to\infty}\psi^{\prime}(\alpha)e^{-\alpha
t}N^{i}_{t},\quad a.s,$ (4.16)

and, let $\mathcal{E}\left(\Xi\right)$ be the random variable defined by

$\mathcal{E}\left(\Xi\right):=\int_{[0,\infty]}\mathcal{E}_{(\xi_{u})}e^{-\alpha
u}\mathds{1}_{\Xi>u}\ \xi(du).$ (4.17)

######  Lemma 4.8 ($L^{2}$ convergence in the general case).

Consider the general frequency spectrum $\left(A(k,t,\Xi)\right)_{k\geq 1}$,
then, for all $k$, $\psi^{\prime}(\alpha)e^{-\alpha t}A(k,t,\Xi)$ converge to
$\mathcal{E}\left(\Xi\right)$ (see 4.17) in $L^{2}$ as $t$ goes to infinity
and

$\lim\limits_{t\to\infty}e^{-\alpha
t}\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t,\Xi)-e^{\alpha
t}\mathcal{E}(\Xi)c_{k}\right)\left(\psi^{\prime}(\alpha)A(l,t,\Xi)-e^{\alpha
t}\mathcal{E}(\Xi)c_{k}\right)\right]=\frac{\alpha}{b}a_{k,l}\int_{\mathbb{R}_{+}}e^{-\alpha
u}\mathbb{P}\left(\Xi>u\right)bdu,$

where the convergence is uniform w.r.t. the random variable $\Xi$. In the case
where $\Xi$ is distributed as $O^{(\beta t)}_{2}$, for $0<\beta<\frac{1}{2}$,
we get

$\lim\limits_{t\to\infty}e^{-\alpha
t}\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t,O^{(\beta
t)}_{2})-e^{\alpha t}\mathcal{E}(O^{(\beta
t)}_{2})c_{k}\right)\left(\psi^{\prime}(\alpha)A(l,t,O^{\beta
t}_{2})-e^{\alpha t}\mathcal{E}(O^{(\beta
t)}_{2})c_{k}\right)\right]=\psi^{\prime}(\alpha)a_{k,l}.$

######  Lemma 4.9 (First moment).

The first moments are asymptotically bounded, that is, for all $k\geq 1$,

$\mathbb{E}\left(\psi^{\prime}(\alpha)A(k,t)(\Xi)-e^{\alpha
t}c_{k}\mathcal{E}(\Xi)\right)\leq\mathcal{O}(1),$

uniformly with respect to the random variable $\Xi$.

######  Lemma 4.10 (Boundedness in the general case.).

Let $k_{1},k_{2},k_{3}$ three positive integers, then

$\mathbb{E}\left[\left|\prod_{i=1}^{3}\frac{\left(\psi^{\prime}(\alpha)A(k_{i},t)-e^{\alpha
t}\mathcal{E}c_{k_{i}}\right)}{e^{\frac{\alpha}{2}t}}\right|\right]=\mathcal{O}\left(1\right),$

uniformly with respect to the random variable $\Xi$.

We do not detail the proofs of these results since they are direct adaptations
of the proofs of Lemmas 6.5, 6.6, and 6.7 of [13].

###  4.2 Proof of the result

The following result is based on the fact that, in the clonal sub-critical
case, the lifetime of a family is expected to be small. It follows that one
can expect that all the family of size $k$ live in different subtrees as soon
as $t>>u$. This is the point of the following lemma.

######  Lemma 4.11.

Suppose that $\alpha<\theta$. If we denote by $\Gamma_{u,t}$ the event,

$\Gamma_{u,t}=\left\\{\text{"there is no family in the population at time t
which is older than }u\text{"}\right\\},$

then, for all $\beta$ in $(0,1-\frac{\alpha}{\theta})$, we have

$\lim\limits_{t\to\infty}\mathbb{P}_{\beta t}\left(\Gamma_{\beta
t,t}\right)=1.$

###### Proof.

The proof of this Lemma, as the calculation of the moments of $A(k,t)$ relies
on the representation of the genealogy of the living population at time $t$ as
a coalescent point process [5]. Moreover, we denote by
$\widetilde{N}^{(t)}_{u}$ the number of living individuals at time $u$ who
have alive descent at time $t$. In [5], we showed that, under
$\mathbb{P}_{t}$, $\widetilde{N}^{(t)}_{u}$ is geometrically distributed with
parameter $\frac{W(t-u)}{W(t)}$.

Now, $\mathds{1}_{\Gamma_{u,t}}$ can be rewritten as

$\mathds{1}_{\Gamma_{u,t}}=\prod_{i=1}^{\widetilde{N}^{(t)}_{u}}\mathds{1}_{\left\\{Z^{i}_{0}(t-u)=0\right\\}},$

where $Z_{0}^{i}(t-u)$ denotes the number of individuals alive at time $t$
descending from the $i$th individual alive at time $u$ and carrying its type
(the clonal type of the sub-CPP). Moreover, from Proposition 4.3 of [5], we
know that that under $\mathbb{P}_{t}$, the family $Z^{(i)}_{0}(t-u)$ is an
i.i.d. family of random variables distributed as $Z_{0}(t-u)$ under
$\mathbb{P}_{t-u}$, and $\widetilde{N}^{(t)}_{u}$ is independent of
$Z^{(i)}_{0}(t-u)$ (still under $\mathbb{P}_{t}$).

Then,

$\mathbb{P}_{t}\left(\Gamma_{t,u}\right)=\mathbb{E}_{t}\left[\mathbb{P}_{t-u}\left(Z_{0}(t-u)=0\right)^{\widetilde{N}^{(t)}_{u}}\right]=\frac{\mathbb{P}_{t-u}\left(Z_{0}(t-u)=0\right)\frac{W(t-u)}{W(t)}}{1-\mathbb{P}_{t-u}\left(Z_{0}(t-u)=0\right)\left(1-\frac{W(t-u)}{W(t)}\right)}.$

Using (2.8), some calculus leads to,

$\mathbb{P}_{t}\left(\Gamma_{t,u}\right)=1-\frac{1}{1+\frac{W_{\theta}(t-u)}{e^{-\theta(t-u)}W(t)}\left(1-\frac{e^{-\theta(t-u)}W(t-u)}{W_{\theta}(t-u)}\right)}.$

Now, since,

$\mathbb{P}_{t}\left(\Gamma_{t,u}\right)=\mathbb{P}_{u}\left(\Gamma_{t,u}\right)\frac{\mathbb{P}\left(N_{u}>0\right)}{\mathbb{P}\left(N_{t}>0\right)}+\frac{\mathbb{P}\left(\Gamma_{t,u},N_{t}=0,N_{u}>0\right)}{\mathbb{P}\left(N_{t}>0\right)},$

taking $u=\beta t$, we obtain, using Lemma 2.2 and

$\mathbb{P}\left(N_{t}=0,N_{\beta t}>0\right)=\mathbb{P}\left(N_{\beta
t}>0\right)-\mathbb{P}\left(N_{t}>0\right)\underset{t\to\infty}{\to}0,$

the desired result. ∎

###### Proof of Theorem 3.1.

Fix $0<u<t$. Note that the event $\Gamma_{u,t}$ of Lemma 4.11 can be rewritten
as

$\mathds{1}_{\Gamma_{u,t}}=\prod_{i=1}^{N_{u}}\mathds{1}_{\left\\{Z^{i}_{0}(t-u,O_{i})=0\right\\}},$
(4.18)

where $Z_{0}^{i}(t-u,O_{i})$ denote the number of individuals alive at time
$t$ carrying the same type as the $i$th alive individual at time $u$, that is
the ancestral family of the splitting constructed from the residual lifetime
of the $i$th individual (see Section 4 in [13]).

Let $K$ be a multi-integer, we denote by $\mathcal{L}^{(K)}$ (resp. $A(K,t)$)
the random vector $\left(\mathcal{L}^{k_{1}},\dots,\mathcal{L}^{k_{N}}\right)$
(resp. $\left(A(k_{1},t),\dots,A(k_{N},t)\right)$) with

$\mathcal{L}_{t}^{k_{i}}=\frac{\psi^{\prime}(\alpha)A(k,t)-c_{k}e^{\alpha
t}\mathcal{E}}{e^{\frac{\alpha}{2}t}}.$

On the event $\Gamma_{u,t}$, we have a.s.,

$A(k_{l},t)=\sum_{i=1}^{N_{u}}A^{(i)}(k_{l},t-u,O_{i}),\quad\forall
l=1,\dots,N,$

where the family $\left(A^{(i)}\left(k_{l},t-u,O_{i}\right)\right)_{i\geq 1}$
stand for the frequency spectrum for each subtree, which are independent from
Lemma 4.1. Hence, using Lemma 4.2,

$\mathcal{L}_{t}^{k_{l}}=\sum_{i=1}^{N_{u}}\frac{\psi^{\prime}(\alpha)A^{(i)}(k_{l},t-u,O_{i})-e^{\alpha(t-u)}\mathcal{E}_{i}(O_{i})c_{k_{l}}}{e^{\frac{\alpha}{2}u}e^{\frac{\alpha}{2}(t-u)}}.$

By Lemma 4.1, that the family $\left(A^{i}(k_{l},t-u,O_{i})\right)_{2\leq
i\leq N_{u}}$ is i.i.d. under $\mathbb{P}_{u}$.

In the sequel, we denote, for all $l$ and $i\geq 1$,

$\tilde{A}^{(i)}\left(k_{l},t-u,O_{i}\right)=\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{l},t-u,O_{i}\right)-e^{\alpha(t-u)}\mathcal{E}_{i}(O_{i})c_{k_{l}}}{e^{\frac{\alpha}{2}(t-u)}}.$

Now, let

$\displaystyle\varphi_{K}\left(\xi\right)$
$\displaystyle:=\mathbb{E}\left[\exp\left(i<\tilde{A}\left(K,t-u,O_{2}\right),\xi>\right)\mathds{1}_{Z^{2}_{0}(t-u,O_{2})=0}\right],$
$\displaystyle\tilde{\varphi}_{K}\left(\xi\right)$
$\displaystyle:=\mathbb{E}\left[\exp\left(i<\tilde{A}\left(K,t-u,O_{1}\right),\xi>\right)\mathds{1}_{Z^{1}_{0}(t-u,O_{1})=0}\right].$

From this point, following closely the proof of Theorem 3.2 of [13]. Taking
$u=\beta$ in $\left(0,\frac{1}{2}\wedge(1-\frac{\alpha}{\theta})\right)$, the
only difficulty is to handle the indicator function
$\mathds{1}_{Z_{0}(t-u,O_{i})>0}$ in the Taylor development of $\varphi_{K}$.
We show how it can be done for one of the second order terms, and leave the
rest of the details to the reader.

It follows from H lder’s inequality that

$\mathbb{E}\left[\left(\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{l},(1-\beta)t,O_{i}\right)-e^{\alpha((1-\beta)t)}\mathcal{E}_{i}(O_{i})c_{k_{l}}}{e^{\frac{\alpha}{2}((1-\beta)t)}}\right)^{2}\mathds{1}_{Z^{2}_{0}((1-\beta)t,O_{2})>0}\right]\\\
\leq\mathbb{E}\left[\left(\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{l},(1-\beta)t,O_{i}\right)-e^{\alpha(1-\beta)t}\mathcal{E}_{i}(O_{i})c_{k_{l}}}{e^{\frac{\alpha}{2}(1-\beta)t}}\right)^{3}\right]^{\frac{2}{3}}\mathbb{P}\left(Z^{2}_{0}((1-\beta)t,O_{2})>0\right)^{\frac{1}{3}},$
(4.19)

from which it follows, using Lemma 4.10, that the r.h.s. of this last
inequality is
$\mathcal{O}\left(\mathbb{P}\left(Z^{2}_{0}(t-u,O_{2})>0\right)^{\frac{1}{3}}\right)$.
Now, using (4.18) and Lemma 4.11, it is easily seen that

$\lim\limits_{t\to\infty}\mathbb{P}\left(Z^{2}_{0}((1-\beta)t,O_{2})>0\right)=0.$

Finally, using Lemma 4.5, we get

$\lim\limits_{t\to\infty}\mathbb{E}\left[\left(\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{l},t-u,O_{i}\right)-e^{\alpha(t-u)}\mathcal{E}_{i}(O_{i})c_{k_{l}}}{e^{\frac{\alpha}{2}(t-u)}}\right)^{2}\mathds{1}_{Z^{2}_{0}(t-u,O_{2})=0}\right]=\psi^{\prime}(\alpha)a_{k,k}.$

These allow us to conclude that

$\lim\limits_{t\to\infty}\mathbb{E}_{\beta
t}\left[e^{i<\mathcal{L}^{(K)}_{t},\xi>}\mathds{1}_{\Gamma_{t}}\right]=\frac{1}{1+\sum_{i,j=1}^{N}\mathcal{M}_{i,j}\
\xi_{i}\xi_{j}},$

where $K_{i,j}$ is given by

$\mathcal{M}_{i,j}:=\psi^{\prime}(\alpha)a_{K_{i},K_{j}},$

with $K$ is the multi-integer $(k_{1},\dots,k_{N})$, and the $a_{l,k}$s are
defined in Lemma 4.5.

To end the proof, note that,

$\left|\mathbb{E}_{\infty}\left[e^{i<\mathcal{L}^{(K)}_{t},\xi>}\right]-\mathbb{E}_{\beta
t}\left[e^{i<\mathcal{L}^{(K)}_{t},\xi>}\mathds{1}_{\Gamma_{\beta
t,t}}\right]\right|\leq\mathbb{E}\left[\left|\frac{\mathds{1}_{\text{NonEx}}}{\mathbb{P}\left(\text{NonEx}\right)}-\frac{\mathds{1}_{N_{\beta
t}>0}\mathds{1}_{\Gamma_{\beta t,t}}}{\mathbb{P}\left(N_{\beta
t}>0\right)}\right|\right]\underset{t\to\infty}{\to}0,$

thanks to Lemma 4.11. ∎

##  5 Proof of Theorem 3.5

Since all the ideas of the proof of this theorem have been developed the
preceding sections, we do not detail all the proof. The only step which needs
clarification is the computation of the covariance matrix of the Laplace limit
law $\mathcal{M}$. According to the proof of Theorem 3.1, it is given by

$\displaystyle\mathcal{M}_{i,j}:=$
$\displaystyle\lim\limits_{t\to\infty}\frac{W(\beta t)}{e^{\alpha\beta
t}}\mathbb{E}\Bigg{[}\left(\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{i},(1-\beta)t,O_{i}\right)-\psi^{\prime}(\alpha)c_{k_{i}}N_{(1-\beta)t}}{e^{\frac{\alpha}{2}((1-\beta)t)}}\right)$
$\displaystyle\times\left(\frac{\psi^{\prime}(\alpha)A^{(i)}\left(k_{j},(1-\beta)t,O_{i}\right)-c_{k_{j}}N_{(1-\beta)t}}{e^{\frac{\alpha}{2}((1-\beta)t)}}\right)\mathds{1}_{Z^{2}_{0}((1-\beta)t,O_{2})>0}\Bigg{]},$

which is equal, thanks to (4.19) and an easy adaptation of Lemma 6.6 in [13],
to

$\displaystyle\mathcal{M}_{i,j}=$
$\displaystyle\lim\limits_{t\to\infty}\frac{b\psi^{\prime}(\alpha)}{\alpha}\frac{W(\beta
t)}{e^{\alpha\beta t}}e^{\alpha t}\mathbb{E}\left[\left(e^{-\alpha
t}A(k_{i},t)-c_{k_{i}}e^{-\alpha t}N_{t}\right)\left(e^{-\alpha
t}A(k_{j},t)-c_{k_{j}}e^{-\alpha t}N_{t}\right)\right].$

So it remains to get the limit of

$e^{\alpha t}\mathbb{E}\left[\left(e^{-\alpha
t}\psi^{\prime}(\alpha)A(k,t)-\psi^{\prime}(\alpha)c_{k}e^{-\alpha
t}N_{t}\right)\left(e^{-\alpha t}\psi^{\prime}(\alpha)A(l,t)-c_{l}e^{-\alpha
t}\psi^{\prime}(\alpha)N_{t}\right)\right],$

as $t$ goes to infinity. We recall that using the calculus made in the proof
of Theorem 6.3 of [5], we have

$\mathbb{E}_{t}A(k,t)N_{t}=2W(t)^{2}c_{k}(t)-2W(t)\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)da+W(t)\int_{[0,t]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]da.$ (5.1)

Moreover, (4.15) entails

$\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}A(k,t)A(l,t)=2W(t)^{2}c_{k}(t)c_{l}(t)+RW(t)+o(e^{-\alpha
t}),$

with

$\displaystyle R:=$
$\displaystyle-\psi^{\prime}(\alpha)\int_{0}^{\infty}2\theta
W(a)^{-1}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[A(l,a)\right]da$
$\displaystyle\quad\quad\quad+\psi^{\prime}(\alpha)\int_{0}^{\infty}2\theta
W(a)^{-1}\mathbb{P}_{a}\left(Z_{0}(a)=l\right)\mathbb{E}_{a}\left[A(k,a)\right]da$
$\displaystyle\quad\quad\quad+\psi^{\prime}(\alpha)\int_{0}^{\infty}\theta
W(a)^{-1}\left(\mathbb{E}_{t}\left[A(k,t)\mathds{1}_{Z_{0}(a)=l}\right]+\mathbb{E}_{t}\left[A(l,t)\mathds{1}_{Z_{0}(a)=k}\right]\right)da.$

These identities allow us to obtain

$\displaystyle\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)\left(A(l,t)-c_{l}N_{t}\right)\right]=2W(t)^{2}c_{k}(t)c_{l}(t)+e^{-\alpha
t}R+o(e^{-\alpha t}),$
$\displaystyle-2c_{l}c_{k}(t)W(t)^{2}+2c_{l}W(t)\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)da-
c_{l}W(t)\int_{[0,t]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle-2c_{k}c_{l}(t)W(t)^{2}+2c_{l}W(t)\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=l\right)da-
c_{k}W(t)\int_{[0,t]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=l}\right]da$
$\displaystyle+c_{k}c_{l}W(t)^{2}\left(2-\frac{1}{W(t)}\right)$
$\displaystyle=2W(t)^{2}\left(c_{k}(t)-c_{l}\right)\left(c_{l}(t)-c_{k}\right)+e^{-\alpha
t}\frac{R}{\psi^{\prime}(\alpha)}+o(e^{-\alpha t}),$
$\displaystyle+2c_{l}W(t)\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)da-
c_{l}W(t)\int_{[0,t]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle+2c_{l}W(t)\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=l\right)da-
c_{k}W(t)\int_{[0,t]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=l}\right]da$
$\displaystyle-c_{k}c_{l}W(t).$

Taking the limit as $t$ goes to infinity leads to

$\displaystyle M_{k,l}:=\lim\limits_{t\to\infty}$
$\displaystyle\psi^{\prime}(\alpha)^{2}e^{-\alpha
t}\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)\left(A(l,t)-c_{l}N_{t}\right)\right]=R$
$\displaystyle+2\psi^{\prime}(\alpha)c_{l}\int_{[0,\infty]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)da-\psi^{\prime}(\alpha)c_{l}\int_{[0,\infty]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle+2\psi^{\prime}(\alpha)c_{l}\int_{[0,\infty]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=l\right)da-\psi^{\prime}(\alpha)c_{k}\int_{[0,\infty]}\theta
W(a)^{-1}\mathbb{E}_{a}\left[N_{a}\mathds{1}_{Z_{0}(a)=l}\right]da$
$\displaystyle-\psi^{\prime}(\alpha)c_{k}c_{l}.$ (5.2)

Finally, since $\mathbb{P}\left(N_{t}>0\right)\sim\frac{\alpha}{b}$,

$\mathcal{M}_{i,j}=M_{k_{i},k_{j}}.$

##  6 Markovian cases

We can get more information on the unknown covariance matrix $K$ in the case
where the life duration distribution is exponential. Our study also cover the
case $\mathbb{P}_{V}=\delta_{\infty}$ (Yule case), although it does not fit
the conditions required by the Theorem 3.1. The reason comes from our method
of calculation for $\mathbb{E}\left[A(k,t)\mathcal{E}\right]$. Let us consider
the filtration $\left(\mathcal{F}_{t}\right)_{t\in\mathbb{R}_{+}}$, where
$\mathcal{F}_{t}$ is the $\sigma$-field generated by the tree truncated above
$t$ and the restriction of the mutation measure on $[0,t)$.

Then $N_{t}$ is Markovian with respect to $\mathcal{F}_{t}$ and for all
positive real numbers $t\leq s$,

$\displaystyle\mathbb{E}\left[A(k,t)N_{s}\mid\mathcal{F}_{t}\right]=A(k,t)N_{t}\mathbb{E}\left[N_{s-t}\right].$

So that,

$\displaystyle\mathbb{E}\left[A(k,t)N_{s}\right]=\mathbb{E}\left[A(k,t)N_{t}\right]\left(W(s-t)-\mathbb{P}_{V}\star
W(s-t)\right).$

By making a renormalization by $e^{-\alpha s}$ and taking the limit as $s$
goes to infinity, we get,

$\mathbb{E}\left[A(k,t)\mathcal{E}\right]=\psi^{\prime}(\alpha)e^{-\alpha
t}\mathbb{E}\left[A(k,t)N_{t}\right],$

since, in the Markovian case, it is known from [8] that

$\frac{\alpha}{b}=\psi^{\prime}(\alpha).$

Suppose first that $d>0$. It follows that,

$\displaystyle\mathbb{E}$
$\displaystyle\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\left(\psi^{\prime}(\alpha)A(l,t)-e^{\alpha
t}c_{l}\mathcal{E}\right)\right]=\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}\left[A(k,t)A(l,t)\right]\mathbb{P}\left(N_{t}>0\right)$
$\displaystyle-
c_{k}\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}\left[A(l,t)N_{t}\right]\mathbb{P}\left(N_{t}>0\right)-c_{l}\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}\left[A(k,t)N_{t}\right]\mathbb{P}\left(N_{t}>0\right)$
$\displaystyle+2\psi^{\prime}(\alpha)e^{2\alpha t}c_{k}c_{l}$

By (4.4),

$\mathbb{P}\left(N_{t}>0\right)=\psi^{\prime}(\alpha)+\psi^{\prime}(\alpha)^{2}\mu
e^{-\alpha t}+o(e^{-\alpha t}),$

so

$\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\left(\psi^{\prime}(\alpha)A(l,t)-e^{\alpha
t}c_{l}\mathcal{E}\right)\right]\\\
=\mathbb{P}\left(N_{t}>0\right)\psi^{\prime}(\alpha)^{2}\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)\left(A(l,t)-c_{l}N_{t}\right)\right]+c_{k}c_{l}\psi^{\prime}(\alpha)\left(2e^{2\alpha
t}-\psi^{\prime}(\alpha)\mathbb{E}_{t}\left[N_{t}^{2}\right]\mathbb{P}\left(N_{t}>0\right)\right).$

Finally, since, using Proposition 2.3,

$\lim\limits_{t\to\infty}e^{-\alpha t}\left(2e^{2\alpha
t}-\psi^{\prime}(\alpha)\mathbb{E}_{t}\left[N_{t}^{2}\right]\mathbb{P}\left(N_{t}>0\right)\right)=\psi^{\prime}(\alpha)\left(1-6\mu\right),$

it follows from (5.2),

$\displaystyle\lim\limits_{t\to\infty}\mathbb{E}$
$\displaystyle\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\left(\psi^{\prime}(\alpha)A(l,t)-e^{\alpha
t}c_{l}\mathcal{E}\right)\right]$
$\displaystyle=\psi^{\prime}(\alpha)M_{k,l}+c_{k}c_{l}\psi^{\prime}(\alpha)^{2}\left(1-6\mu\right)=\psi^{\prime}(\alpha)M_{k,l}+c_{k}c_{l}\psi^{\prime}(\alpha)^{2}\left(1-6\frac{d}{\alpha}\right),$

using that $\mu=\frac{1}{b\mathbb{E}V-1}$. In the Yule case, an easy
adaptation of the preceding proof leads to

$\displaystyle\lim\limits_{t\to\infty}\mathbb{E}\left[\left(\psi^{\prime}(\alpha)A(k,t)-e^{\alpha
t}c_{k}\mathcal{E}\right)\left(\psi^{\prime}(\alpha)A(l,t)-e^{\alpha
t}c_{l}\mathcal{E}\right)\right]=M_{k,l}+c_{k}c_{l}.$

##  7 Numerical studies

The purpose of this section is to analyze our approximation method and the
estimation of the error by virtue of numerical experiments. There are several
practical difficulties appearing when one tries to perform such study.

The first problem, which involves no conceptual difficulties, lies only on the
implementation of the formulas appearing in Theorems 3.1, 3.3 and 3.5. In
particular, the computation of the moments of type
$\mathbb{E}[A(k,t)\mathds{1}_{Z_{0}(t)=l}]$ are particularly complicated (see
Proposition 5.4 in [5]).

Another difficulty is to obtain numerical approximations of the scale
functions $W$ and $W_{\theta}$. For instance, these functions appear in the
computation of the covariance matrix of Theorems 3.3 and 3.5 or when one wants
to simulate the _coalescent point process_. To obtain such approximations, we
need to apply numerically the Laplace inversion operator to the functions
$\frac{1}{\psi}$ and $\frac{1}{\psi_{\theta}}$.

Unfortunately, the Laplace numerical inversion is a rather difficult problem
(see for instance [1] or [2]) which is often computationally expensive. As a
consequence, the computational cost of performing multiple numerical
integration involving $W$ or $W_{\theta}$ can be important when done with a
crude method. Moreover, these methods presents rough numerical instabilities
when the original function is exponentially increasing (inverting
$\lambda\to\frac{1}{1-\lambda}$, whose inverse is $x\to e^{x}$, is already a
tough numerical problem).

For all these reasons, we provide with this work a Matlab toolbox which handle
all these difficulties and allows users to perform numerical experiments
without having to take care of these issues.

In this whole section, we are interested in the approximation of the frequency
spectrum at a fixed time $t$ by the sequence $N_{t}(c_{k})_{k\geq 1}$ (we
recall that $c_{k}$ was defined in equation (2.13)). As a consequence, the
error in this approximation are computed thanks to Theorem 3.5. The parameters
of the model are set as follows:

  * •

$\mathbb{P}_{V}$ is a Rice distribution with shape parameter 1 and scale
parameter 1.

  * •

$b=1$.

  * •

$\theta=1$.

For such parameters $\alpha$ approximately equals to $0.5$. Figure 3 shows the
evolution of the frequency spectrum (for $k$ between $1$ and $10$) through
time. The different quantities seem to growth exponentially with rate $\alpha$
with a time-shift which depend on $k$. An interesting open question would be
to understand the behavior of these shifts.

Figure 3: A simulation of the evolution of the frequency spectrum under the
given model.

In order to stress our methods of approximation, the first idea is to look to
the renormalized frequency spectrum $\left(\frac{A(k,t)}{c_{k}}\right)_{k\geq
1}$ which is expected to look like $(N_{t},\ t\in\mathbb{R}_{+})$. As showed
in Figure 4, the approximation seems to be quite accurate for $k=1,2$.
However, a more quantitative analysis is required. Figure 5 shows the absolute
error in the approximation of $A(1,t)$ by $c_{1}N_{t}$. This error is a little
disappointing since it since to diverge when $t$ goes to infinity. However,
even if, according to Figure 5, the absolute error at time $20$ if of order
$10^{3}$, the relative error shows that this error is quite small with respect
to the value of $A(1,20)$.

Figure 4: Evolution of the renormalized frequency spectrum
$(A(k,t)/c_{k})_{k\geq 1}$ under the given model.

Figure 5: Absolute error (left picture) and relative (right picture) in the
approximation of $A(1,t)$ by $c_{1}N_{t}$.

Another question is about the speed of convergence in the central limit
theorem stated in Theorem 3.5. The red curve of Figure 6 shows the density of
the Laplace distribution given in Theorem 3.5 in the case of $A(1,t)$ whereas
the blue histogram shows the distribution of
$\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ for $t=10$
($\alpha t\sim 5$ and $\mathbb{E}_{t}[N_{t}]\sim 300$) from 10000 simulations.
This Figure highlights the fact that even if the taken time $t$ is quite small
the distribution
$\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ seems
already close to the limiting distribution. Figure 7 shows the same kind of
behavior in the multidimensional case. To be more quantitative, Figure 8 shows
the evolution in time of the distance between the density of limit
distribution given in Theorem 3.5 and a kernel estimation of the distribution
of $\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ (the
estimation is made from $10000$ simulations at each time). This suggest an
exponential rate of convergence in Theorem 3.5. In the view of Figure 8, one
may think that Berry-Essen type results for Theorem 3.5 would be quite
interesting, in particular to understand how the speed of convergence is
related to choice of the parameters.

Figure 6: Distribution of the renormalized error and expected limit
distribution given by our CLT. Figure 7: Joint distribution of the
renormalized error (left figure) and expected limit distribution (right
figure) given by our CLT. Figure 8: Estimation of the rate of convergence in
$L^{2}$ norm.

Another interesting question which could be partially probed by simulation is
the study of the behavior of the error in the clonal supercritical case.
Figure 9 shows a kernel estimation (from $10000$ simulation) of the density of
$\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ in the
clonal supercritical case ($\theta=0.2<\alpha$). Figure 9 suggest a totally
different behavior with a limit distribution which is asymmetric with respect
to $0$. In particular, in the view of the shape of the distribution, one could
conjecture that the limit is a skew stable distribution.

Figure 9: Kernel estimate of the probability density function of the limit
distribution in the clonal supercritical case.

To end this section, let us goes back to one of the motivation of this work.
The following discussion dot not claim to be rigorous and is essentially
formal. We recall that the Extended Haplotype Homozygosity (EHH) can be used
to detect positive selection in a population [18]. In particular, the behavior
of the frequency spectrum in this model gives a standard for the behavior of a
subpopulation carrying a common allele under neutral evolution. In order to
have a rigorous model to describe this phenomenon, we need to introduce a new
mutation measure which is different from the one given in Section 2. We define
it directly on the CPP but this could be equivalently defined on the splitting
tree. So let $\mathcal{P}$ be a Poisson random measure on
$[0,t]\times\mathbb{N}\times\mathbb{R}_{+}$ with intensity measure
$\lambda\otimes C\otimes\lambda$, where $C$ is the counting measure on
$\mathbb{N}$, then, for any mutation rate $\theta$ in $\mathbb{R}_{+}$, we
define the $\theta$-mutation random measure $\mathcal{N}_{\theta}$ by

$\mathcal{N}_{\theta}\left(A\times B\right)=\int_{A\times
B\times[0,\theta]}\mathds{1}_{H_{i}>t-a}\mathds{1}_{i<\mathcal{N}_{t}}\mathcal{P}\left(di,da,dx\right),$

where, as before, an atom at $(a,i)$ means that the $i$th branch experiences a
mutation at time $t-a$. This construction allows to increase the mutation in
consistent manner. This allows to model the type of an individual at a
distance $x$ (such that the mutation rate is a function of $x$) from the core
haplotype (we refer the reader to [18] for more details). Now, following [5],
we can define the frequency spectrum at mutation rate $\theta$ by

$A^{\theta}(k,t)=\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z_{0}(i,a)=k}\mathcal{N}_{\theta}(di,da),$

where $Z_{0}(i,a)$ is the number of individual at time $t$ carrying the type
of the $i$th individual at time $t-a$ (see [5] for more details). Let us also
define $Z_{0}^{\theta}(t)$ the number of individuals carrying the type of the
first individual at time $0$ when the mutation measure is given by
$\mathcal{N}_{\theta}$. As expected, the allelic partition of the population
becomes thinner as $\theta$ growth.

Now, the definition of the EHH $G_{\theta}(t)$ is the probability that two
uniformly sampled individuals in the population have the same type, that is

$G_{t}(\theta)=\frac{Z^{\theta}_{0}(t)(Z^{\theta}_{0}(t)-1)+\sum_{k\geq
1}k(k-1)A^{\theta}(k,t)}{N_{t}(N_{t}-1)}.$

Using that

$N_{t}=Z^{\theta}_{0}(t)+\sum_{k\geq 1}kA^{\theta}(k,t),$

this rewrite

$G_{t}(\theta)=\frac{(N_{t}-\sum_{k\geq 1}kA^{\theta}(k,t))(N_{t}-\sum_{k\geq
1}kA^{\theta}(k,t)-1)+\sum_{k\geq 1}k(k-1)A^{\theta}(k,t)}{N_{t}(N_{t}-1)}.$

Finally, using the approximation

$(A(k,t))_{k\geq 1}\approx(c_{k})_{k\geq 1}N_{t}$

proposed in this work, one could expect that

$G_{t}(\theta)\approx\frac{\sum_{k\geq
1}k(k-1)c_{k}}{N_{t}}=\frac{\int_{0}^{\infty}2\theta e^{-\theta
x}(W_{\theta}(x)-1)dx}{N_{t}}.$

We stress the fact that the above expression make sens only in the clonal
subcritical case (in the other cases the integral in note finite). Now, we can
look at the accuracy of this approximation in view of numerical simulation.
Figure 10 shows the value of the EHH (when $\theta$ increase) from a
simulation of the model (blue curve) and the one obtained using our
approximation (red curve). In view of Figure 10, the approximation seems
pretty accurate. In order to be more quantitative, Figure 11 shows the
relative error between the EHH and its approximation for one simulation. This
shows that the error, as least for sufficiently large $\theta$, remains under
$8$%. To end, let us highlight that Theorem 3.5 can be used to give confidence
intervals for fixed $\theta$ but in order the construct tests of selection
from curves like these of Figure 10 one would need to have functional CLT in
long time for the process $((A^{\theta}(k,t)-c^{\theta}_{k}N_{t})_{k\geq 1},\
\theta\in\mathbb{R}_{+})$.

Figure 10: Extended Haplotype Homozygosity (EHH) with the given approximation
(red curve) and from simulated data (blue curve) Figure 11: Relative error in
the approximation of the EHH

##  Appendix A A bit of renewal theory

The purpose of this part is to recall some facts on renewal equations borrowed
from [10]. Let $h:\mathbb{R}\to\mathbb{R}$ be a function bounded on finite
intervals with support in $\mathbb{R}_{+}$ and $\Gamma$ a probability measure
on $\mathbb{R}_{+}$. The equation

$F(t)=\int_{\mathbb{R}_{+}}F(t-s)\Gamma(ds)+h(t),$

called a renewal equation, is known to admit a unique solution finite on
bounded interval.

Here, our interest is focused on the asymptotic behavior of $F$. We said that
the function $h$ is DRI (directly Riemann integrable) if for any $\delta>0$,
the quantities

$\delta\sum_{i=0}^{n}\sup_{t\in[\delta i,\delta(i+1))}f(t)$

and

$\delta\sum_{i=0}^{n}\inf_{t\in[\delta i,\delta(i+1))}f(t)$

converge as $n$ goes to infinity respectively to some real number
$I_{sup}^{\delta}$ and $I_{inf}^{\delta}$, and

$\lim\limits_{\delta\to 0}I_{sup}^{\delta}=\lim\limits_{\delta\to
0}I_{inf}^{\delta}<\infty.$

In the sequel, we use the two following criteria for the DRI property:

######  Lemma A.1.

Let $h$ a function as defined previously. If $h$ satisfies one of the next two
conditions, then $h$ is DRI:

  * 1.

$h$ is non-negative decreasing and classically Riemann integrable on
$\mathbb{R}_{+}$,

  * 2.

$h$ is càdlàg and bounded by a DRI function.

We can now state the next result, which is constantly used in the sequel.

######  Theorem A.2.

Suppose that $\Gamma$ is non-lattice, and $h$ is DRI, then

$\lim\limits_{t\to\infty}F(t)=\gamma\int_{\mathbb{R}_{+}}h(s)ds,$

with

$\gamma:=\left(\int_{\mathbb{R}_{+}}s\ \Gamma(ds)\right)^{-1},$

if the above integral is finite, and zero otherwise.

######  Remark A.3.

In particular, if we suppose that $\Gamma$ is a measure with mass lower than
$1$, and that there exists a constant $\alpha\geq 0$ such that

$\int_{\mathbb{R}_{+}}e^{\alpha t}\Gamma(dt)=1,$

then, one can perform the change a measure

$\widetilde{\Gamma}(dt)=e^{\alpha t}\Gamma(dt),$

in order to apply Theorem A.2 to a new renewal equation to obtain the
asymptotic behavior of $F$. (See [10] for details). This method is also used
in the sequel.

##  Appendix B Formula for the fourth moment of the error

######  Lemma B.1.

$\displaystyle\mathbb{E}_{t}$
$\displaystyle\left[\left(A(k,t)-c_{k}N_{t}\right)^{4}\right]=4\int_{[0,t]}\theta\frac{W(t)}{W(a)}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)^{3}\right]da$
$\displaystyle+48\int_{[0,t]}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}N_{a}A(k,a)\right]\mathbb{E}_{a}\left[\left(c_{k}N_{a}-A(k,a)\right)\right]da$
$\displaystyle+24\int_{[0,t]}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}N_{a}^{2}\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]da$
$\displaystyle+24\int_{[0,t]}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)^{2}\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]da$
$\displaystyle+8\int_{[0,t]}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)^{3}\right]da$
$\displaystyle+48\int_{[0,t]}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)^{2}\right]da$
$\displaystyle+72\int_{[0,t]}\theta\frac{W(t)^{3}}{W(a)^{3}}\left(1-\frac{W(a)}{W(t)}\right)^{2}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{2}da$
$\displaystyle+72\int_{[0,t]}\theta\frac{W(t)^{3}}{W(a)^{3}}\left(1-\frac{W(a)}{W(t)}\right)^{2}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)^{2}\right]\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle+96\int_{[0,t]}\theta\frac{W(t)^{4}}{W(a)^{4}}\left(1-\frac{W(a)}{W(t)}\right)^{3}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{3}da+c_{k}^{4}\mathbb{E}_{t}N_{t}^{4}$

###### Proof.

The proof of this Lemma lies on the calculation of the expectation of each
term in the development of

$\left(A(k,t)-c_{k}N_{t}\right)^{4}.$

To make this, we intensively use the relation (2.10) and the method developed
in [5]. We begin by computing

$\mathbb{E}_{t}\left[A(k,t)^{4}\right].$

Formula (2.10) gives us,

$\displaystyle A(k,t)^{4}=$ $\displaystyle
4\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k_{i}}\sum_{u_{1:3}=1}^{N^{(t)}_{t-a}}\prod_{\underset{i\neq
j}{j=1}}^{3}A^{(u_{j})}(k,a)\mathcal{N}\left(da,di\right)$ $\displaystyle=$
$\displaystyle
4\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}A^{i}(k,a)A^{i}(k,a)A^{i}(k,a)\mathcal{N}(da,di)$
$\displaystyle+4\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}\sum_{\underset{j_{1}\neq
j_{2}\neq j_{3}\neq
i}{j_{1},j_{2},j_{3}=1}}^{N^{(t)}_{t-a}}A^{j_{1}}(k,a)A^{j_{2}}(k,a)A^{j_{3}}(k,a)\mathcal{N}(da,di)$
$\displaystyle+12\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}A^{i}(k,a)A^{i}(k,a)\sum_{j=1,j\neq
i}^{N^{(t)}_{t-a}}A^{j}(k,a)\mathcal{N}(da,di)$
$\displaystyle+4\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}\sum_{j=1,j\neq
i}^{N^{(t)}_{t-a}}A^{j}(k,a)^{3}\mathcal{N}(da,di)$
$\displaystyle+12\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}A^{i}(k,a)\sum_{j_{1},j_{2}=1,j_{1}\neq
j_{2}\neq i}^{N^{(t)}_{t-a}}A^{j_{1}}(k,a)A^{j_{2}}(k,a)\mathcal{N}(da,di)$
$\displaystyle+24\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}A^{i}(k,a)\sum_{j_{1}=1,j_{1}\neq
i}^{N^{(t)}_{t-a}}A^{j_{1}}(k,a)A^{j_{1}}(k,a)\mathcal{N}(da,di)$
$\displaystyle+12\int_{[0,t]\times\mathbb{N}}\mathds{1}_{Z^{i}_{0}(a)=k}\sum_{j_{1},j_{2}=1,j_{1}\neq
j_{2}\neq
i}^{N^{(t)}_{t-a}}A^{j_{1}}(k,a)^{2}A^{j_{2}}(k,a)\mathcal{N}(da,di).$ (B.1)

The decomposition of the sum in form

$\sum_{u_{1:3}=1}^{N^{(t)}_{t-a}},$

has then been made to distinguish independence properties in our calculation.
Actually, as soon as, $i\neq j$, $A^{i}(k,a)$ is independent from $A^{i}(k,a)$
(see [5] for details). It is essential to note that the expectation of these
integrals with respect to the random measure $\mathcal{N}$ are all calculated
thanks to Theorem 3.1 of [5]. So, taking the expectation now leads to,

$\displaystyle\mathbb{E}_{t}\left[A(k,t)^{4}\right]=$ $\displaystyle
4\int_{[0,t]}\theta\mathbb{E}_{a}\left[N^{(t)}_{t-a}\right]\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)^{3}\right]\theta
da$
$\displaystyle+4\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(4)}\right]\mathbb{E}_{a}\left[A(k,a)\right]^{3}da$
$\displaystyle+12\int_{[0,t]}\theta\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)^{2}\right]\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(2)}\right]\mathbb{E}_{a}\left[A(k,a)\right]da$
$\displaystyle+4\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(2)}\right]\mathbb{E}_{a}\left[A(k,a)^{3}\right]da$
$\displaystyle+12\int_{[0,t]}\theta\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)\right]\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(3)}\right]\mathbb{E}_{a}\left[A(k,a)\right]^{2}da$
$\displaystyle+24\int_{[0,t]}\theta\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}A(k,a)\right]\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(2)}\right]\mathbb{E}_{a}\left[A(k,a)^{2}\right]da$
$\displaystyle+12\int_{[0,t]}\theta\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(N^{(t)}_{t-a}\right)_{(3)}\right]\mathbb{E}_{a}\left[A(k,a)^{2}\right]\mathbb{E}_{a}\left[A(k,a)\right]da.$

Using the same method for all the other terms and that, for any positive real
number $a$ lower than $t$,

$N_{t}=\sum_{i=1}^{N^{(t)}_{t-a}}N^{(i)}_{a},$

we get Lemma B.1 by reassembling similar terms together. The last term is
obtained using the geometric distribution of $N_{t}$ under $\mathbb{P}_{t}$. ∎

##  Appendix C Boundedness of the fourth moment

######  Lemma C.1.

We begin the proof of the boundedness of the fourth moment by some estimates.

$\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)\right]=\mathcal{O}\left(e^{-(\theta-\alpha)t}\right),$
(i)
$\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)^{3}\right]=\mathcal{O}\left(W(t)^{2}\right),$
(ii)
$\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)^{2}\right]=\mathcal{O}\left(W(t)\right),$
(iii) $\mathbb{E}_{t}N_{t}^{n}=\mathcal{O}(e^{n\alpha t}),\quad
n\in\mathbb{N}^{*},$ (iv)
$\mathbb{P}_{t}\left(Z_{0}(t)=k\right)=\mathcal{O}(e^{(\alpha-\theta)t}).$ (v)

###### Proof.

Relation (i) is easily obtained using the expectation of $N_{t}$ and $A(k,t)$
using (2.11), (2.13) and the behaviour of $W$ provided by Proposition 2.3. The
relation (iii) has been obtained in the proof of Theorem 6.1 in [5]. The two
last relations are easily obtained from (2.4), (2.8) and Lemma 2.2. The
relation (ii) is obtained using the following estimation,

$\left|\mathbb{E}_{t}\left[\left(A(k,a)-c_{k}N_{a}\right)^{3}\right]\right|\leq\mathbb{E}_{t}\left[N_{a}\left(A(k,a)-c_{k}N_{a}\right)^{2}\right].$

We begin the proof by computing the r.h.s. of the previous inequality using
the same techniques as in Appendix A.

$\displaystyle\mathbb{E}\left[A(k,t)^{2}N_{t}\right]$
$\displaystyle=2\int_{0}^{t}\theta\frac{W(t)}{W(a)}\mathbb{E}\left[N_{a}A(k,a)\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]\mathbb{E}\left[A(k,a)\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}\left[A(k,a)\mathds{1}_{Z_{0}(a)=k}\right]\mathbb{E}\left[N_{a}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[A(k,a)N_{a}\right]da$
$\displaystyle+$ $\displaystyle
12\int_{0}^{t}\theta\frac{W(t)^{3}}{W(a)^{3}}\left(1-\frac{W(a)}{W(t)}\right)^{2}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[A(k,a)\right]\mathbb{E}\left[N_{a}\right]da.$

$\displaystyle 2\mathbb{E}\left[A(k,t)N_{t}^{2}\right]$
$\displaystyle=2\int_{0}^{t}\theta\frac{W(t)}{W(a)}\mathbb{E}\left[N_{a}^{2}\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle+$ $\displaystyle
8\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]\mathbb{E}\left[N_{a}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[N_{a}^{2}\right]da$
$\displaystyle+$ $\displaystyle
12\int_{0}^{t}\theta\frac{W(t)^{3}}{W(a)^{3}}\left(1-\frac{W(a)}{W(t)}\right)^{2}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[N_{a}\right]^{2}da.$

Finally,

$\displaystyle\mathbb{E}\left[N_{t}\left(A(k,t)-c_{k}N_{t}\right)^{2}\right]$
$\displaystyle=2\int_{0}^{t}\theta\frac{W(t)}{W(a)}\mathbb{E}\left[N_{a}\left(A(k,a)-c_{k}N_{a}\right)\mathds{1}_{Z_{0}(a)=k}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}\left[N_{a}\mathds{1}_{Z_{0}(a)=k}\right]\mathbb{E}\left[A(k,a)-c_{k}N_{a}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{E}\left[\left(A(k,a)-c_{k}N_{a}\right)\mathds{1}_{Z_{0}(a)=k}\right]\mathbb{E}\left[N_{a}\right]da$
$\displaystyle+$ $\displaystyle
4\int_{0}^{t}\theta\frac{W(t)^{2}}{W(a)^{2}}\left(1-\frac{W(a)}{W(t)}\right)\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[N_{a}\left(A(k,a)-c_{k}N_{a}\right)\right]da$
$\displaystyle+$ $\displaystyle
12\int_{0}^{t}\theta\frac{W(t)^{3}}{W(a)^{3}}\left(1-\frac{W(a)}{W(t)}\right)^{2}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}\left[N_{a}\right]\mathbb{E}\left[A(k,a)-c_{k}N_{a}\right]da$
$\displaystyle+$ $\displaystyle c_{k}^{2}\mathbb{E}_{t}N_{t}^{3}.$

Now, an analysis similar to the one of Lemma 4.7 leads to the result. ∎

###### Proof of Lemma 4.7.

The ideas of the proof, is to analyses one to one every terms of the
expression of

$\mathbb{E}_{t}\left[\left(A(k,t)-c_{k}N_{t}\right)^{4}\right],$

given by Lemma B.1 using Lemma C.1 to show that they behave as
$\mathcal{O}\left(W(t)^{2}\right)$. Since the ideas are the same for every
terms, we just give a few examples.

First of all, we consider

$\int_{[0,t]}\frac{W(t)}{W(a)}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)^{3}\right]da.$

Using Lemma C.1 (ii), we have

$\int_{[0,t]}\frac{W(t)}{W(a)}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)^{3}\right]da=\mathcal{O}\left(W(t)^{2}\right).\\\
$

Now take the term

$\int_{[0,t]}\frac{W(t)^{2}}{W(a)^{2}}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}N_{a}^{2}\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]da,$

we have from Lemma C.1 (i) and (iv),

$\int_{[0,t]}\frac{W(t)^{2}}{W(a)^{2}}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}N_{a}^{2}\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]da\leq\int_{[0,t]}\frac{W(t)^{2}}{W(a)^{2}}\mathbb{E}_{a}\left[N_{a}^{2}\right]e^{-(\theta-\alpha)a}da=\mathcal{O}\left(W(t)^{2}\right).\\\
$

Every term in $W(t)$ or $W(t)^{2}$ are treated this way. Now, we consider the
term in $W(t)^{4}$ which is

$I:=96\int_{[0,t]}\frac{W(t)^{4}}{W(a)^{4}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{3}da+24W(t)^{4}c_{k}^{4},$

since $N_{t}$ is geometrically distributed under $\mathbb{P}_{t}$, and that

$\mathbb{E}_{t}N_{t}^{4}=24W(t)^{4}-36W(t)^{3}+\mathcal{O}(W(t)^{2}).$ (C.1)

On the other hand, using the law of $Z_{0}(t)$ given by (2.8) and the
expectation of $A(k,t)$ given by (2.11) (under $\mathbb{P}_{t}$), we have,

$\displaystyle 96$
$\displaystyle\int_{[0,t]}\frac{W(t)^{4}}{W(a)^{4}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{3}da$
$\displaystyle=-96W(t)^{4}\int_{0}^{t}\frac{\theta e^{-\theta
a}}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\left(\int_{0}^{a}\frac{\theta
e^{-\theta s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}\
ds\right)^{3}\ da$ $\displaystyle=-24W(t)^{4}\left(\int_{0}^{t}\frac{\theta
e^{-\theta a}}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\
da\right)^{4}.$

Finally,

$I=24W(t)^{4}\left(\int_{t}^{\infty}\frac{\theta e^{-\theta
a}}{W_{\theta}(a)^{2}}\left(1-\frac{1}{W_{\theta}(a)}\right)^{k-1}\
da\right)^{4}=\mathcal{O}\left(W(t)^{4}e^{-4\theta t}\right)=o(1).$

The last example is the most technical and relies with the term in $W(t)^{3}$,
which is, using (C.1) and Lemma B.1,

$\displaystyle J:=$ $\displaystyle
72\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)\right]\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{2}da$
$\displaystyle+72\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)^{2}\right]\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle-288\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{3}da-36c_{k}^{4}W(t)^{3}.$

On the other hand, using the calculus made in the proof of Theorem 6.3 of [5],
we have

$\displaystyle\mathbb{E}_{a}$
$\displaystyle\left[\left(A(k,a)-c_{k}N_{a}\right)^{2}\right]$
$\displaystyle=$ $\displaystyle
4\int_{[0,a]}\frac{W(a)^{2}}{W(s)^{2}}\left(1-\frac{W(s)}{W(a)}\right)\mathbb{P}_{s}\left(Z_{0}(s)=k\right)\mathbb{E}_{a}\left(A(k,s)-c_{k}N_{s}\right)ds$
$\displaystyle+2\int_{[0,a]}\frac{W(s)}{W(a)}\
\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(s)=k}\left(A(k,s)-c_{k}N_{s}\right)\right]ds+c_{k}^{2}W(a)^{2}\left(2-\frac{1}{W(a)}\right).$

Substituting this last expression in $J$ leads to

$\displaystyle J$
$\displaystyle=-144\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)\right]\int_{[a,\infty]}\frac{\mathbb{P}\left(Z_{0}(a)=k\right)}{W(s)^{2}}\mathbb{E}_{a}\left[\left(A(k,s)-c_{k}N_{s}\right)\right]dsda$
$\displaystyle+144W(t)^{3}\int_{[0,t]}\frac{1}{W(a)}\
\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)\right]\int_{[a,t]}\frac{1}{W(s)^{2}}\mathbb{P}_{s}\left(Z_{0}(s)=k\right)\mathbb{E}_{a}\left[A(k,s)-N_{s}c_{k}\right]da$
$\displaystyle-144c_{k}^{2}\int_{[0,t]}\frac{W(t)^{3}}{W(a)}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle+144\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{P}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]^{3}da$
$\displaystyle-288\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{2}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\int_{[0,a]}\frac{1}{W(s)}\mathbb{P}_{s}\left(Z_{0}(s)=k\right)\mathbb{E}_{a}\left(A(k,s)-c_{k}N_{s}\right)ds\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle+72\int_{[0,t]}\frac{W(t)^{3}}{W(a)}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)c_{k}^{2}\left(2-\frac{1}{W(a)}\right)\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle-288\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{P}_{a}\left(Z_{0}(a)=k\right)\mathbb{E}_{a}\left[\left(A(k,a)-c_{k}N_{a}\right)\right]^{3}da-36c_{k}^{4}W(t)^{3}.$

Using many times that,

$\displaystyle\int_{[0,t]}$
$\displaystyle\frac{\theta\mathbb{P}\left(Z_{0}(a)=k\right)}{W(s)^{2}}\mathbb{E}_{a}\left[\left(A(k,s)-c_{k}N_{s}\right)\right]ds$
$\displaystyle=$ $\displaystyle-\int_{[0,t]}\frac{\theta e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}\int_{[s,\infty]}\frac{\theta
e^{-\theta
u}}{W_{\theta}(u)^{2}}\left(1-\frac{1}{W_{\theta}(u)}\right)^{k-1}duds$
$\displaystyle=$
$\displaystyle\frac{c_{k}^{2}}{2}-\frac{1}{2}\left(\int_{[t,\infty]}\frac{\theta
e^{-\theta
s}}{W_{\theta}(s)^{2}}\left(1-\frac{1}{W_{\theta}(s)}\right)^{k-1}ds\right)^{2},$

thanks to (2.8), (2.11), and (2.6), we finally get

$\displaystyle J=$
$\displaystyle-144\left(c_{k}^{2}-c_{k}(t)^{2}\right)\int_{[0,t]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{E}_{a}\left[\mathds{1}_{Z_{0}(a)=k}\left(A(k,a)-c_{k}N_{a}\right)\right]da$
$\displaystyle+36W(t)^{3}\left(c_{k}^{2}\left(\int_{[t,\infty]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]^{3}da\right)^{2}-\left(\int_{[t,\infty]}\frac{W(t)^{3}}{W(a)^{3}}\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]^{3}da\right)^{4}\right)$
$\displaystyle+144\left(c_{k}-c_{k}(t)\right)^{2}\int_{[0,t]}\frac{W(t)^{3}}{W(a)}\mathbb{E}_{a}\left[A(k,a)-N_{a}c_{k}\right]da$
$\displaystyle+36W(t)^{3}\left(c_{k}-c_{k}(t)\right)^{4}.$

This shows that $J$ is $\mathcal{O}\left(W(t)^{2}\right)$. ∎

## References

  * [1] Joseph Abate, Gagan L Choudhury, and Ward Whitt.  An introduction to numerical transform inversion and its application to probability models.  In Computational probability, pages 257–323. Springer, 2000.
  * [2] Joseph Abate and Ward Whitt.  A unified framework for numerically inverting laplace transforms.  INFORMS Journal on Computing, 18(4):408–421, 2006.
  * [3] K. B. Athreya and P. E. Ney.  Branching processes.  Dover Publications, Inc., Mineola, NY, 2004.  Reprint of the 1972 original [Springer, New York; MR0373040].
  * [4] Jean Bertoin.  The structure of the allelic partition of the total population for Galton-Watson processes with neutral mutations.  Ann. Probab., 37(4):1502–1523, 2009.
  * [5] Nicolas Champagnat and Henry Benoit.  Moments of the frequency spectrum of a splitting tree with neutral poissonian mutations.  Electron. J. Probab., 21:34 pp., 2016.
  * [6] Nicolas Champagnat and Amaury Lambert.  Splitting trees with neutral Poissonian mutations I: Small families.  Stochastic Process. Appl., 122(3):1003–1033, 2012.
  * [7] Nicolas Champagnat and Amaury Lambert.  Splitting trees with neutral Poissonian mutations II: Largest and oldest families.  Stochastic Process. Appl., 123(4):1368–1414, 2013.
  * [8] Nicolas Champagnat, Amaury Lambert, and Mathieu Richard.  Birth and death processes with neutral mutations.  Int. J. Stoch. Anal., pages Art. ID 569081, 20, 2012.
  * [9] Warren J. Ewens.  Mathematical population genetics. I, volume 27 of Interdisciplinary Applied Mathematics.  Springer-Verlag, New York, second edition, 2004.  Theoretical introduction.
  * [10] William Feller.  An introduction to probability theory and its applications. Vol. II. Second edition. John Wiley & Sons, Inc., New York-London-Sydney, 1971\.
  * [11] J. Geiger and G. Kersting.  Depth-first search of random trees, and Poisson point processes.  In Classical and modern branching processes (Minneapolis, MN, 1994), volume 84 of IMA Vol. Math. Appl., pages 111–126. Springer, New York, 1997.
  * [12] R. C. Griffiths and Anthony G. Pakes.  An infinite-alleles version of the simple branching process.  Adv. in Appl. Probab., 20(3):489–524, 1988.
  * [13] B. Henry.  Central limit theorem for supercritical binary homogeneous Crump-Mode-Jagers processes.  to appear in ESAIM:Probability and Statistics, November 2016.
  * [14] Amaury Lambert.  The contour of splitting trees is a Lévy process.  Ann. Probab., 38(1):348–395, 2010.
  * [15] Amaury Lambert, Lea Popovic, et al.  The coalescent point process of branching trees.  The Annals of Applied Probability, 23(1):99–144, 2013.
  * [16] Amaury Lambert and Pieter Trapman.  Splitting trees stopped when the first clock rings and vervaat’s transformation.  Journal of Applied Probability, 50(01):208–227, 2013.
  * [17] Mathieu Richard.  Arbres, Processus de branchement non Markoviens et processus de Lévy.  Thèse de doctorat, Université Pierre et Marie Curie, Paris 6.
  * [18] Pardis C Sabeti, David E Reich, John M Higgins, Haninah ZP Levine, Daniel J Richter, Stephen F Schaffner, Stacey B Gabriel, Jill V Platko, Nick J Patterson, Gavin J McDonald, et al.  Detecting recent positive selection in the human genome from haplotype structure.  Nature, 419(6909):832–837, 2002.