# Ergodic behaviour of a multi-type growth-fragmentation process modelling the mycelial network of a filamentous fungus

Milica Tomasevic, Vincent Bansaye, Amandine Véber

► To cite this version:

Milica Tomasevic, Vincent Bansaye, Amandine Véber. Ergodic behaviour of a multi-type growth-fragmentation process modelling the mycelial network of a filamentous fungus. ESAIM: Probability and Statistics, 2022, 26, pp.397-435. 10.1051/ps/2022013. hal-03087196v3

HAL Id: hal-03087196

https://hal.science/hal-03087196v3

Submitted on 20 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ergodic behaviour of a multi-type growth-fragmentation process modelling the mycelial network of a filamentous fungus

Milica Tomašević *

Vincent Bansaye †

Amandine Véber ‡

## Abstract

In this work, we introduce a stochastic growth-fragmentation model for the expansion of the network of filaments, or *mycelium*, of a filamentous fungus. In this model, each individual is described by a discrete type $e \in \{0,1\}$ indicating whether the individual corresponds to an internal or terminal segment of filament, and a continuous trait $x \ge 0$ corresponding to the length of this segment. The length of internal segments cannot grow, while the length of terminal segments increases at a deterministic speed $v$. Both types of individuals/segments branch according to a type-dependent mechanism.

After constructing the stochastic bi-type growth-fragmentation process of interest, we analyse the corresponding mean measure (or first moment semigroup). We show that its ergodic behaviour is, as expected, governed by the maximal eigenelements. In the long run, the total mass of the mean measure increases exponentially fast while the type-dependent density in trait converges to an explicit distribution $N$, independent of the initial condition, at some exponential speed. We then obtain a law of large numbers that relates the long term behaviour of the stochastic process to the limiting distribution $N$. In the particular model we consider, which depends on only 3 parameters, all the quantities needed to describe this asymptotic behaviour are explicit, which paves the way for parameter inference based on data collected in lab experiments.

**Key words and phrases:** *Branching processes, growth-fragmentation systems, limit theorems, spectral analysis.*

**MSC2020 subject classification:** 60J80, 35B40, 60J85, 92-10.

## 1 Introduction

Filamentous fungi are complex expanding organisms that are omnipresent in nature. They form filamentous structures known as *hyphae*. These filaments grow and branch to create potentially huge networks called *mycelia*, sometimes covering up to a few square kilometers. To feed the whole mycelium, hyphae tamper with their environment by decomposing the dead organic matter, making its chemical components available to the next generation of organisms (including the fungus itself). Thereby, filamentous fungi play a key role in the

<footnote>* Corresponding author. CMAP, Ecole Polytechnique, CNRS, I.P. Paris, route de Saclay 91128 Palaiseau, France. milica.tomasevic@polytechnique.edu.</footnote>

<footnote>†CMAP, Ecole Polytechnique, CNRS, I.P. Paris, route de Saclay 91128 Palaiseau, France. vincent.bansaye@polytechnique.edu.</footnote>

<footnote>‡MAP5, Université de Paris Cité, CNRS, 45 rue des Saints-Pères 75006 Paris, France. amandine.veber@parisdescartes.fr.</footnote>
functioning of natural ecosystems. They are also able to quickly respond to local threats
such as attacks by predators, physical obstacles, or noxious local conditions, through an
efficient chemical communication along the hyphal network, leading to its partial reorgan-
isation or to the reorientation of the growth capacity of the mycelium away from danger
[10, 25]. All these characteristics made these species appealing to the biochemical (and
in particular, pharmaceutical) industry, in which they are now routinely used to catalyse
various reactions and produce different types of metabolites. See the review paper [8]
for numerous examples and for a careful discussion of the relation between the extent of
mycelial branching and metabolite production.

In this paper, we aim at understanding the basic growth properties of the mycelium
in a given species of filamentous fungi, leaving aside the complex interactions with its
environment and ecosystem mentioned above. More precisely, assuming that the fungus
grows in homogeneous conditions and away from predators or pests (e.g., on a Petri dish
in a lab experiment), we want to identify simple descriptors that characterise the growth
of the fungus and allow us to quantify the impact of various forms of stresses (nutrient
depletion, pH, ... ) on the mycelial growth and structure. These descriptors should be
robust enough to allow the comparison of different fungus species.

The expansion of the hyphal network rests upon several biological processes. The first
one is the growth of “primary” hyphae, in a more or less radial way. These hyphae extend
the area already covered by the mycelium, exploring the environment in search for new
sources of nutrients. They grow in numbers by branching in two at their tips (or *apexes* –
we shall later speak of *apical* branching) at some rate. As in a spider web, these primary
hyphae serve as a backbone for “secondary” hyphae, that branch off from the primary
structure (approximately uniformly along the existing hyphae – we shall later speak of
*lateral* branching). The secondary hyphae increase the density of the network by growing
in different directions and by themselves branching both laterally and at their apexes.
A third process is the fusion of two hyphae when they cross, called *anastomosis*. This
phenomenon improves the connectivity of the mycelium, as it creates shortcuts for the
diffusion of molecules along the network of filaments. Note that hyphae do not necessarily
merge when they cross, even when they are constrained to evolve in two spatial dimensions
as in lab conditions. Alternatively, they may simply bypass each other and keep growing in
different directions; anastomosis represents only a fraction of the outcomes of the crossing
events and these crossings may in fact occur less frequently in nature (in which fungi grow
in three dimensions) than on the two dimensional surface of a Petri dish.

Understanding the basic growth properties and the branching structure of the network
of hyphae that results from them will be the first step before engaging in a more detailed
modelling of the fungal growth taking into account flows of nutrients and chemical signalling
along the network, which will be the object of future work. These questions have already
been the object of a lot of attention, and the current state of the art in models of mycelial
growth is substantial. A significant part of it relies on graph theory and consists in the
statistical analysis of experimental quantitative data (number of internal nodes, of apexes,
hyphal length, ..., of mycelia grown in laboratory) [8, 20, 26]. The major difficulties in these
approaches, on which progress is still being made, is to set up a high-quality recording of
this multi-scale growth dynamics (starting from a spore of a few micrometres and ending
when the mycelium covers the few square centimetres of a Petri dish), and to be able to
extract the topological network information of interest thanks to semi- or fully automated
post-processing tools able to resolve most of the ambiguities present in the images (note
that the number of branch points observable in these images can reach 10⁵, rendering
node identification “with the eye” clearly unfeasible). We refer to [20] for more details on current challenges in this area. More mechanistic approaches gave rise to a variety of spatially explicit stochastic models for the spatial spread of fungal mycelia, either lattice-based or lattice-free, in which hyphae grow in length by colonising neighbouring (free) locations, branch at some rate to give birth to a new hypha colonising free locations in another direction, and sometimes merge with another hypha when they cross. See [11] for a review of such models and [22] and references therein for more recent work. Despite their clear mathematical formulation, to our knowledge no analytical results on the long-term growth properties of these stochastic networks have been derived, in particular due to the difficulty of handling spatial interactions such as anastomosis and self-avoidance. The analysis instead relies on intensive simulations, exploring the space of parameters to find families of parameters producing patterns that match the observations, or looking for phase transitions in the mycelial growth pattern.

On a much finer scale, most models zoom in on the tip of a single hypha to understand the mechanisms triggering its extension (see, e.g., [4, 33]). In contrast, other models focus on much larger spatial and temporal scales and describe the interaction between filamentous fungi and their environments in natural conditions, mostly in an aggregated way : the whole mycelium is assimilated with a single scalar quantity, its biomass, and a system of ordinary differential equations describes the circulation of nutrients and chemical molecules between the fungus and the environment and its effect on the growth and degradation of the fungus biomass and on different characteristics of the environment (see, e.g., [28]). More recently, spatially explicit models were introduced, based on reaction-diffusion partial differential equations [23] or based on a system of stochastic differential equations (encoding the behaviour of each hypha, including anastomosis and self-avoidance) and its mean-field deterministic limit [16]. These models are in the same vein as models of tumor-induced angiogenesis (see, e.g., [15]) and allow in particular to study global quantities such as the stationary shape and speed of the invasion front formed by the whole mycelium on the spatial scale of observation with the naked eye.

In this work, we instead focus on an intermediate spatial scale, observable in lab experiments thanks to the previously mentioned technology for the automated recording and analysis of movies such as the one developed at LIED (University of Paris), which was used in [20] to measure several quantities (exponential growth rates of the number of free ends of filaments – “open” ends – and of the number of internal nodes in the network – “closed” ends – in particular) analogous to those which will naturally appear in our analysis below. Because we shall neglect anastomosis in order to keep a tractable model as a first exploration tool in what follows, some care will be needed when doing parameter inference based on the model developed below. This issue will be discussed in Section 6.

Let us now describe our approach. It is based on two strong assumptions which are mostly motivated by our aim to understand the exponential growth behaviour observed in [20] through a simple but informative half-mechanistic, half-statistical model :

(i) In completely homogeneous conditions and over small space- and time-scales, the spatial organisation of the hyphae does not (really) matter. That is, since the mycelium naturally spreads over the available space, we shall make the approximation that every piece of filament evolves in the same conditions as the others. Here we neglect the depletion of food due to the high density of hyphae around the origin (or centre) of the mycelium, which is a reasonable choice when considering short timescales;

(ii) Anastomosis (the fusion of crossing hyphae) does not need to be explicitly modelled
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Schematic representation of a network of filaments. Black dots represent branch points (except the left-most dot which corresponds to the start point of the network), while open ends that keep on elongating are depicted by triangles. In this example, the network is made of 15 segments of filaments, of which 7 are internal (i.e., lie between two branch points) and 8 are terminal (i.e., lie between a branch point and the open end of a filament).</figcaption></figure>

and its impact can be incorporated via an appropriate statistical treatment of the data when doing parameter inference. See the section on *Apex, node and length growth dynamics* in [20] and Section 6 below for more details on this point.

The second assumption may look like a surprising modelling choice, but it has the paramount advantage that distinct (pieces of) filaments will not interact with each other in the model, allowing us to encode the mycelium as a *branching process*.

Each individual in our branching process corresponds to a segment of filament lying between two branching points (internal, or *closed*, segment), or between a branching point and the extremity of the filament (terminal, or *open*, segment). See Figure 1. An individual is represented by a pair $(\mathfrak{e}, x)$, where $x \ge 0$ is the current length of the corresponding segment and $\mathfrak{e} \in \{0, 1\}$ encodes whether the segment is open ($\mathfrak{e} = 1$) or closed ($\mathfrak{e} = 0$). More precisely, the space in which the pairs $(\mathfrak{e}, x)$ take their values is

$$S := (\{0\} \times (0, \infty)) \cup (\{1\} \times [0, \infty)), \qquad (1.1)$$

excluding closed segments of length 0 for mathematical convenience. Indeed, in the dynamics defined below the state $(0, 0)$ is an absorbing state and individuals with these characteristics will a.s. never be produced by the fragmentation of other individuals (see Remark 2.1). Therefore, should we include $(0, 0)$ in the definition of $S$, the mass of such individuals in the population would remain constant equal to its initial value, and taking this mass into account would only create an artificial particular case to be considered in each step of our analysis.

At every time $t \ge 0$, the set of all segments constituting the network is fully described by the following point measure on $S$

$$Z_t := \sum_{u \in V_t} \delta_{(\mathfrak{e}^u, x_t^u)}, \qquad (1.2)$$

where $V_t$ denotes the indexing set of the individuals alive at time $t$ and $z_t^u = (\mathfrak{e}^u, x_t^u)$ denotes the characteristics at time $t$ of individual $u \in V_t$. Note that this representation of the population of segments at any given time does not allow us to infer who is hooked up with whom in the network, but it will be sufficient for our purposes (cf. Assumption (i)).
<figure><img src="image_3.png" /><figcaption>**Figure 2:** (a) Apical branching. An open segment branches into two new filaments at its end. The open end of the segment therefore closes (and the segment becomes “closed” itself) and two open segments, initially of length 0, are created and start growing in length. (b) Lateral branching. A new open segment of initial length zero branches off from an existing segment, which can be open or closed, at a location which is uniformly distributed along the current length of the existing segment. This branching event fragments the existing length into 2 segments (one necessarily closed, and one of the same type as the fragmented segment), and adds a third segment which is open and initially of length zero.</figcaption></figure>

We write $\mathcal{M}_p(S)$ for the space of all finite point measures on $S$ and we endow it with the topology of weak convergence. We shall also use the standard notation, for $\nu = \sum_{i=1}^n \delta_{z_i}$ and $\varphi$ a measurable function on $S$,

$$ \langle \nu, \varphi \rangle := \int_S \varphi(z) \nu(dz) = \sum_{i=1}^n \varphi(z_i). $$

Let us fix $v, b_1, b_2 \in (\mathbb{R}_+^*)^3$. The dynamics of the process $(\mathcal{Z}_t)_{t \ge 0}$ are as follows:

(a) **Elongation.** Open segments elongate deterministically at speed $v$, while closed segments cannot grow. More precisely, for every $s, t \ge 0$ and $u \in V_s$, conditionally on individual $u$ not being involved in a branching event during the time interval $[s, s+t)$, we have for all $r \in [0, t)$,

$$ x_{s+r}^u = x_s^u + \mathcal{E}^u vr. \qquad (1.3) $$

(b) **Apical branching.** Each open segment branches “at its apex” at rate $b_1$. That is, every extant individual $u$ such that $\mathcal{E}^u = 1$, independently of each other and at rate $b_1$, is removed from the population and replaced by three new individuals: one closed individual of the same length as the “parent”, therefore with characteristics $(0, x_{t-}^u)$ (where $t$ is the time of the branching event), and two open individuals of length 0 (and thus both with characteristics $(1, 0)$). See Figure 2(a).

(c) **Lateral branching.** Every segment (open or closed) of length $x > 0$ branches laterally at rate $b_2x$ and the branch point is chosen uniformly at random along the segment. That is, for every $t \ge 0$ and every $u \in V_{t-}$ such that $x_{t-}^u > 0$, at the instantaneous rate $b_2x_{t-}^u$ individual $u$ is removed from the population and is replaced by three new individuals: one closed individual of length $\alpha^u x_{t-}^u$, one individual with first characteristics $\mathcal{E}^u$ (i.e., open if the parent was open, or closed if the parent was closed) and length $(1 - \alpha^u)x_{t-}^u$, and finally an open individual of length 0, where $\alpha^u$ is an independent draw from a uniform distribution over $[0, 1]$. See Figure 2(b). All individuals branch “laterally” independently of each other and independently of the apical branching events described in (b).

In the above, removing (resp., adding) an individual naturally translates into removing (resp., adding) the corresponding atom in $\mathcal{Z}_t$ at the time of the branching event. Inspired
by this description, let us introduce the following operator $\mathcal{G}$, on which the martingale problem satisfied by $(\mathcal{Z}_t)_{t \ge 0}$ will be based. Let $C_b^1(\mathbb{R})$ stand for the set of all bounded functions on $\mathbb{R}$ of class $C_b^1$ with bounded derivative, and let $C_b^1(S)$ stand for the set of all bounded continuous functions on $S$ with bounded and continuous derivatives w.r.t. the variable $x$. For every $F \in C_b^1(\mathbb{R})$, $f \in C_b^1(S)$, let the function $F_f$ be defined by

$$F_f(\nu) := F(\langle \nu, f \rangle), \quad \nu \in \mathcal{M}_p(S), \qquad (1.4)$$

and define for all such functions $F_f$ and all $\nu \in \mathcal{M}_p(S)$:

$$\begin{aligned}
\mathcal{G}F_f(\nu) := F'(\langle \nu, f \rangle) & \int_S \mathfrak{e}v \frac{\partial f}{\partial x}(\mathfrak{e}, x)\nu(d\mathfrak{e}, dx) \qquad (1.5) \\
& + b_1 \int_S \mathfrak{e} \{F(\langle \nu, f \rangle - f(\mathfrak{e}, x) + f(0, x) + 2f(1, 0)) - F(\langle \nu, f \rangle)\} \nu(d\mathfrak{e}, dx) \\
& + b_2 \int_S x \int_{(0,1)} \{F(\langle \nu, f \rangle - f(\mathfrak{e}, x) + f(0, (1-\alpha)x) + f(\mathfrak{e}, \alpha x) + f(1, 0)) \\
& \qquad \qquad \qquad - F(\langle \nu, f \rangle)\} d\alpha \nu(d\mathfrak{e}, dx).
\end{aligned}$$

Note that if $\nu$ gives positive mass to the point $(1, 0) \in S$, by convention we use the right limit $f(0, 0+)$ to give a sense to the a priori undefined term $f(0, 0)$ appearing in the second integral on the r.h.s of (1.5). In Remark 2.1, we shall argue that the Lebesgue measure of the set of times $t$ at which the measure describing the current state of the population has an atom at $(1, 0)$ is zero, so that the chosen convention is unimportant.

In Section 2, we follow [29] and construct a process $(\bar{Z}_t)_{t \ge 0}$ on a larger space in which the genealogical relationship between individuals is retained through the standard Ulam-Harris-Neveu encoding $\mathcal{U}$. This is the result of Theorem 2.2. If we then restrict our attention to its marginal over $S$ and write $(\mathcal{Z}_t)_{t \ge 0}$ for the resulting $\mathcal{M}_p(S)$-valued process (see Equation (2.6) for a more precise definition), we have the following property, under a first moment assumption which is enough for our purposes. It is also proved in Section 2. Let $p_l: S \to \mathbb{R}_+$ be the projector on the “length” coordinate, defined by $p_l(\mathfrak{e}, x) = x$ for all $(\mathfrak{e}, x) \in S$.

**Proposition 1.1.** Let $Z^0$ be a random variable with values in $\mathcal{M}_p(S)$ such that

$$\mathbb{E}[\langle Z^0, 1 \rangle] < \infty \quad \text{and} \quad \mathbb{E}[\langle Z^0, p_l \rangle] < \infty. \qquad (1.6)$$

Let $(\mathcal{Z}_t)_{t \ge 0}$ be the càdlàg $\mathcal{M}_p(S)$-valued process constructed in (2.5)-(2.6), starting at $\mathcal{Z}_0 = Z^0$. Then for every $F \in C_b^1(\mathbb{R})$ and $f \in C_b^1(S)$,

$$\left(F_f(\mathcal{Z}_t) - F_f(\mathcal{Z}_0) - \int_0^t \mathcal{G}F_f(\mathcal{Z}_s)ds\right)_{t \ge 0} \qquad (1.7)$$

is a martingale.

To be more precise, we have tacitly assumed that $\bar{Z}^0$ is constructed from $Z^0$ by giving a label 1, 2, ..., $\langle Z^0, 1 \rangle$ to the atoms of $Z^0$ to obtain a point measure on $\mathcal{U} \times S$. See Section 2 for more details. Observe that Theorem 2.2 ensures that $\mathcal{Z}_t$ is well-defined at any time $t \ge 0$, which is not obvious from the informal description of the process in terms of growth and branching events. Indeed, the total branching rate at time $t$ is proportional to the number of open individuals $\langle \mathcal{Z}_t, \mathbb{1}_{\{\mathfrak{e}=1\}} \rangle$ and to the total length $\langle \mathcal{Z}_t, p_l \rangle$. Since the
total length process itself increases between the branching times at a speed proportional to the number of open individuals, we need to check that no explosion occurs. We also derive useful bounds on the expectation of the total number of individuals and of the total length at time t in Lemma 2.3. From now on, we work with the process $(Z_t)_{t \ge 0}$ of Proposition 1.1.

What we have just defined is a stochastic growth-fragmentation model, in which the growth term simply corresponds to the elongation of each open individual at a fixed linear speed while the branching of an individual corresponds to its fragmentation “into three pieces” (some of length 0) in such a way that the total length is conserved. Growth-fragmentation models are now common, in particular in the literature on partial differential equations where they have been extensively studied. See [32] for a classic and [5, 13, 21] for more recent examples covering different biological phenomena. Other approaches based on probabilistic representations of the solutions to growth-fragmentation equations appeared more recently, see in particular [7, 9] and references therein for the non-conservative case of which our model is an example. Potential applications include the evolution of age-structured populations (where age is the continuous individual trait that grows linearly and is “fragmented” into one individual with the same age and one individual with age 0 at the event of a birth), the growth of bacterial populations (in which the length, or another continuous individual trait of a bacterium, grows during its lifetime and is split between the two offspring bacteria resulting from a division event), as well as more general fragmentation phenomena (e.g., of polymers) in which potentially more than two individuals come out of a fragmentation event and the sharing of the “parental” trait between the “offspring” may not be conservative (with the production of dust, for instance). Of particular interest is the long-term behaviour of the population size and trait distribution. A huge literature is devoted to this question, and we refer to the introduction of [30] for a comprehensive overview. In general, what is shown is that there exists an exponent $\lambda \in \mathbb{R}$ (called the *Malthusian exponent*) and a stationary profile $N$ on the trait space such that the density $n_t(x)$ of individuals of trait $x$ at time $t$ behaves like

$$n_t(x) \approx e^{\lambda t} \langle n_0, \psi \rangle N(x) \quad \text{as } t \to \infty, \qquad (1.8)$$

where $\psi$ is a function characterising the impact of the initial condition. The approximation (1.8) is made rigorous by considering the appropriate function space and by proving the convergence of $e^{-\lambda t} n_t$ to $\langle n_0, \psi \rangle N$ in this space. Although multidimensional continuous traits (age and size, ...) are sometimes considered, the case of additional non-evolving discrete types that influence the individual growth and branching properties (like our type $\mathcal{C} \in \{0, 1\}$) is less common and we could only find a few studies dealing with quiescent and proliferative cells in models for tumour growth (see, e.g., Section 6 in [1], or [12]).

Such deterministic approaches apply when one considers very large populations, already distributed on the trait space according to some continuous density at the origin of time, for which we suspect that the stochasticity inherent to the branching or fragmentation dynamics only plays a minor role and the average behaviour of the population is sufficient to understand how the trait distribution evolves in time. To deal with initially small populations, or to justify the growth-fragmentation equation at the population level as being the large-population limit of a model at the microscopic level of individual dynamics, individual-based stochastic models are particularly appropriate tools. The classical connection between the branching process $(Z_t)_{t \ge 0}$ and deterministic growth-fragmentation equations is made through the *mean measure* (or first moment semigroup) defined as follows. Let $\mathcal{M}_f(S)$ be the space of all finite measures on $S$ (also endowed with the topology of weak convergence). For every $t \ge 0$, let $n_t \in \mathcal{M}_f(S)$ be defined by: for every bounded
measurable function $f$ on $S$,

$$ \langle n_t, f \rangle := \mathbb{E}[\langle \mathcal{Z}_t, f \rangle]. \qquad (1.9) $$

Note that these quantities are well-defined for all $t$ thanks to the bound on $\mathbb{E}[\langle \mathcal{Z}_t, 1 \rangle]$ obtained in Lemma 2.3. Note also that the definition of $n_t$ depends on the distribution of $\mathcal{Z}_0$, although for now we do not report this dependence in the notation for simplicity. To ease the statement of our next results, let us decompose each $n_t$ as follows: for every bounded measurable $f: S \to \mathbb{R}$,

$$ \langle n_t, f \rangle = \int_{\mathbb{R}_+} f(1, x)n_t^1(dx) + \int_{\mathbb{R}_+^*} f(0, x)n_t^0(dx), \qquad (1.10) $$

where the measure $n_t^1$ on $\mathbb{R}_+$ (resp., $n_t^0$ on $\mathbb{R}_+^*$) is uniquely defined by the set of equations (1.10) written for all $f$ such that $f(\mathfrak{e}, x) = 0$ whenever $\mathfrak{e} = 0$ (resp., $\mathfrak{e} = 1$).

In Section 3, we use the martingale problem formulation (1.7) to show that if $n_0^1$ and $n_0^0$ both admit a density with respect to Lebesgue measure, then this property also holds true at any later time and, furthermore, these densities solve a system of growth-fragmentation equations. More precisely, we prove the following result.

**Proposition 1.2.** Suppose that the assumptions of Proposition 1.1 are satisfied, and that the measures $n_0^1$ and $n_0^0$ defined by (1.10) (with $t = 0$) are absolutely continuous with respect to Lebesgue measure on $\mathbb{R}_+$ and $\mathbb{R}_+^*$, respectively. Then for every $t \ge 0$, $n_t^1$ and $n_t^0$ are also absolutely continuous with respect to Lebesgue measure. Furthermore, abusing notation and writing

$$ n_t^1(dx) = n_t^1(x) dx \quad \text{and} \quad n_t^0(dx) = n_t^0(x) dx, \qquad (1.11) $$

then the couple $(n_t^0, n_t^1)_{t \ge 0}$ is a weak solution to the following system: For all $x > 0$,

$$ \left\{ \begin{array}{l} \frac{\partial}{\partial t} n_t^1(x) + v \frac{\partial}{\partial x} n_t^1(x) + (b_1 + b_2x)n_t^1(x) = b_2 \int_x^\infty n_t^1(y)dy, \\ \\ vn_t^1(0) = 2b_1 \int_0^\infty n_t^1(y)dy + b_2 \int_0^\infty y(n_t^1(y) + n_t^0(y))dy, \\ \\ \frac{\partial}{\partial t} n_t^0(x) + b_2xn_t^0(x) = b_1n_t^1(x) + 2b_2 \int_x^\infty n_t^0(y)dy + b_2 \int_x^\infty n_t^1(y)dy. \end{array} \right. \qquad (1.12) $$

The notion of weak solution that we use here will be made precise in the proof.

Once this system has been derived, we may use it to understand the long-term behaviour of the mean measure (instead of the full stochastic process, to start with). To do so, let us first observe that if we scale time by a factor $1/v$ and consider $(n_{t/v})_{t \ge 0}$, then all the above remains true but elongation now happens at speed $\tilde{v} = 1$, apical branching at rate $\tilde{b}_1 = b_1/v$ and lateral branching at rate $\tilde{b}_2x := (b_2/v)x$. Therefore, to ease the notation and without loss of generality, we now suppose that $v = 1$. Second, let us introduce the following functions, which will be needed for our convergence theorem below. Let $\lambda > 0$ be the unique positive solution to

$$ b_1 + \frac{b_2}{\lambda} = \lambda. \qquad (1.13) $$

We shall see later that $\lambda$ is the maximal eigenvalue of the spectral problem associated to (1.12), and we now turn to the corresponding eigenvector.

Let $N_1: \mathbb{R}_+ \to \mathbb{R}_+$ and $N_0: \mathbb{R}_+^* \to \mathbb{R}_+$ be defined by

$$ N_1(x) = (b_1 + b_2x + \lambda)e^{-\int_0^x (b_1 + b_2y + \lambda)dy}, \qquad (1.14) $$
and

$$N_0(x) = \frac{N_1(x)}{(b_2x + \lambda)^2} (b_2 + b_1(b_2x + \lambda)) + \frac{b_2e^{-\int_0^x (b_1+b_2y+\lambda)dy}}{(b_2x + \lambda)^3} (2b_2 + b_1(b_2x + \lambda)). \quad (1.15)$$

We shall show in Proposition 4.1 that $N_1$ and $N_0$ are probability densities (i.e., they integrate to 1). Let us also define the functions $\psi$ and $V$ on $S$ as follows:

$$\psi(\mathfrak{e}, x) = \frac{\lambda^2}{\lambda^2 + b_2} \left( \mathfrak{e} + \frac{b_2}{\lambda} x \right), \quad V(\mathfrak{e}, x) = \psi(\mathfrak{e}, x) + 1 + x^2. \quad (1.16)$$

These quantities may look fairly mysterious at the moment. We shall see in Section 4.1 that they are the eigenelement of the adjoint problem associated to (1.12) (see Proposition 4.1) and a Lyapunov function useful to control the dynamics starting from large initial values, respectively. We can now formulate the following theorem.

**Theorem 1.3.** There exist $C, w > 0$ such that for any solution $(n_t^0, n_t^1)_{t \ge 0}$ of (1.12) satisfying $\langle n_0, V \rangle < \infty$, we have for all $t \ge 0$

$$\|e^{-\lambda t} n_t^1 - \langle n_0, \psi \rangle N_1\|_{L^1(\mathbb{R}_+)} + \|e^{-\lambda t} n_t^0 - \langle n_0, \psi \rangle N_0\|_{L^1(\mathbb{R}_+^*)} \le C e^{-wt} \langle n_0, V \rangle. \quad (1.17)$$

Theorem 1.3 is in fact a corollary of the finer convergence theorem formulated in Section 4.2, namely Theorem 4.3, in which the convergence is shown to occur in weighted total variation norm in the appropriate functional space. Since its statement requires another round of heavy notation, we only provide here a more reader-friendly version. The proof of Theorem 4.3 uses an analogue of Harris' ergodic theorem for non-conservative semigroups (i.e., whose total mass is not conserved through time) derived from growth-fragmentation equations. More precisely, we use Theorem 2.1 in [7] and most of Section 4.2 will be devoted to proving that the assumptions of this theorem are satisfied by $(n_t)_{t \ge 0}$.

This approach was also used in the recent paper [19] in which individuals (bacteria) can be of two types with different growth parameters: for both types of individuals, the trait considered grows exponentially fast but at two different rates $\alpha_0, \alpha_1 > 0$, and the branching/fragmentation rate is common to all individuals and is trait-dependent. During each fragmentation event, the length of the “parent” is split between the two offspring in fixed proportions $\theta_0, \theta_1 = 1 - \theta_0$. Because of the very quick elongation of both types of individuals, together with the assumption that the (positive) branching rate tends to infinity as the individual length goes to infinity, it is natural (although not at all easy to prove) that the same form of convergence (1.8) as in similar systems with only one type of individuals should occur in their framework, and indeed this constitutes the main results of [19]. In contrast, in our case closed individuals do not elongate and the length of open individuals increases rather slowly (linearly with time). Since the branching rate of closed individuals is proportional to their lengths, it is not at all obvious that the subpopulation of closed individuals will develop quickly enough that the whole population size increases exponentially fast. Furthermore, smaller and smaller closed individuals may accumulate, preventing the length distribution within the population from stabilising. Therefore, the convergence stated in Theorem 1.3 is a somehow more surprising example of the robustness of the asymptotic behaviour of growth-fragmentation equations.

Finally, we can build on Theorem 1.3 to obtain the long-term behaviour of our multi-type growth-fragmentation stochastic model $(Z_t)_{t \ge 0}$ through a law of large numbers. In-
deed, let us define the following measure $N$ on $S$, in the same spirit as the decomposition (1.10): for every bounded measurable $f: S \to \mathbb{R}$,

$$ \langle N, f \rangle = \int_{\mathbb{R}_+} f(1, x) N_1(dx) + \int_{\mathbb{R}_+^*} f(0, x) N_0(dx), $$

where $N_1$ and $N_0$ are the probability distributions defined in (1.14) and (1.15). Our last result is the following.

**Theorem 1.4.** Assume that the conditions stated in (1.6) are satisfied. Then for every measurable function $f: S \to \mathbb{R}$ such that $\sup_{(\mathfrak{c},x) \in S} |f(\mathfrak{c}, x)|/(1+x) < \infty$, the following convergence in probability holds:

$$ \lim_{t \to \infty} \frac{\langle \mathcal{Z}_t, f \rangle}{\langle \mathcal{Z}_t, 1 \rangle} = \frac{1}{2} \langle N, f \rangle. $$

The proof of this law large numbers follows a well-established strategy. It exploits the classical martingale associated to the eigenvalue $\lambda$, $(Y_t)_{t \ge 0} := (\exp(-\lambda t)\langle \mathcal{Z}_t, \psi \rangle)_{t \ge 0}$, and the decorrelation properties of the genealogy of the most recent commun ancestors. More precisely, following [6, 18] and using the spectral gap guaranteed by Theorem 4.3, we shall prove a convergence in $L^2$ via a formula for forks (or many-to-two formula) and simultaneously estimate the speed of convergence. More generally, there is a long story of law of large numbers for multitype branching processes. The $L \log L$ criterion guaranteeing the non-degenerescence of the limit of the martingale $Y$ and the a.s. convergence of the proportions of each type of individuals within the population (without moment assumptions) hold for the case of a finite number of types [3, 27]. These results admit various extensions in infinite dimensions (see, e.g., [2, 24]), involving some additional spectral or moment assumptions. In particular, [2] guarantees a.s. convergence under an additional hypothesis of uniformity of the approximation of the asymptotic profile of the first moment semigroup by eigenelements. Adapting the arguments of [2, 24] would probably allow us to obtain the a.s. convergence of the ratios considered in Theorem 1.4. Our result is weaker, but the approach chosen here for the estimation of the first moment semigroup and the law of large numbers both cover the class of test functions in which we are interested and should be well suited for future relevant extensions of this work. Indeed, the techniques of proof used in this work should easily extend to more generalised forms for the fragmentation and growth rates as well as to varying environments, and should allow us to estimate the speed of convergence of the estimators of interest in such models.

Roughly speaking, Theorem 1.4 means that the empirical distribution $\mathcal{Z}_t/\langle \mathcal{Z}_t, 1 \rangle$ is close to the probability distribution $(1/2)N$ when $t$ is large. In view of our application, if we could consider that the mycelium grew only by branching and elongation (disregarding anastomosis), this would allow us to set up a statistical method to infer the three parameters of the model from experimental data on the type- and length-distribution of the segments of mycelium observed at some large time $t$. In particular, a direct consequence of Theorem 1.4 and of the relation

$$ \lambda \int_0^\infty x(N_1(x) + N_0(x))dx = \int_0^\infty N_1(x)dx $$

that we shall derive at the end of Appendix A.1 is the following. Recall the notation $p_l$ for the projector on the “length” coordinate, and let $p_e: S \to \{0, 1\}$ denote the projector
on the “type” coordinate defined by $p_e(c, x) = c$, for all $(c, x) \in S$. For $t \ge 0$, let us define the following estimator of $\lambda$:

$$ \hat{\Lambda}_t := \frac{\langle \mathcal{Z}_t, p_e \rangle}{\langle \mathcal{Z}_t, p_l \rangle}. \qquad (1.18) $$

Applying Theorem 1.4 to $p_l$ and $p_e$, we obtain that $\hat{\Lambda}_t$ converges in probability to $\lambda$ as $t \to \infty$. Hence, in the absence of anastomosis, the exponential growth rate of the (total) number of segments can be simply estimated by the ratio of the number of “open” segments (or extremal pieces of filaments) to the total length of the network at some large time $t$. In Section 6 we shall briefly discuss how the fusion of filaments distorts these expectations and the directions we shall pursue in future work to get around this issue.

Finally, note that the theorems from [29] and [7] on which our existence and long-term convergence results hinge are very general and would allow us to consider various generalisations of our type-dependent growth and fragmentation mechanisms. However, in the particular case considered here all the quantities and functions of interest are explicit (a nice property which is bound to fail for most generalisations) and in future work we intend to apply this precise model to experimental data, in collaboration with our colleagues from the NEMATIC<fnref n="1" /> research project. We shall therefore stick to this particular model, which should be seen as a simple but characteristic example of what may be done with other multi-type growth-fragmentation dynamics.

The rest of the paper is organised as follows. In Section 2 we construct the $M_p(S)$-valued process of interest, $(\mathcal{Z}_t)_{t \ge 0}$, and provide useful bounds on the expected total size $\langle \mathcal{Z}_t, 1 \rangle$ and total length $\langle \mathcal{Z}_t, p_l \rangle$ of the system at any time $t \ge 0$. In Section 3, we prove Proposition 1.2. In Section 4, we analyse the long time behaviour of the mean measure. First, we state the spectral problem associated to the system (1.12) and its adjoint and provide explicit solutions to both problems (Section 4.1). Second, in Section 4.2, we prove the finer convergence Theorem 4.3, of which Theorem 1.3 is a consequence. Theorem 1.4 is proved in Section 5, and the challenges related to parameter inference are discussed in Section 6. Finally, we prove several technical results in the appendices. In Appendix A.1, we show that the eigenvalues of the above mentionned spectral problem are solutions to Equation (1.13). In Appendix A.2, we check that the explicit values we provided as the solutions to the spectral problem indeed satisfy it.

## 2 Construction of the process $(\mathcal{Z}_t)_{t \ge 0}$

We work on a probability space $(\Omega, \mathcal{F}, \mathbb{P})$ rich enough to accommodate all the objects we need below. Following the method used in [29], we first construct a process $(\bar{\mathcal{Z}}_t)_{t \ge 0}$ keeping track of the different lines of descent. To do so, we use the standard Ulam notation to identify each individual:

$$ \mathcal{U} = \bigcup_{n \in \mathbb{N}} \{\mathbb{N} \times (\{1, 2, 3\})^n\}. \qquad (2.1) $$

For a given initial state of the population $Z_0$ with $I_0$ individuals (or atoms in $Z_0$), we label its atoms by $1, \dots, I_0$. The offspring of an individual $u \in \mathcal{U}$ are denoted by $u1, u2, u3$.

We now proceed as in Section 2 of [29], with a few simplifications due to our particular framework. First, the growth and branching dynamics are homogeneous in time and consequently, in their notation we may take $\mathcal{X} = \mathcal{Y} = S$ (we do not need to add a last

<footnote n="1">¹Growing and branching networks: Analysis, modelling and simulation of multi-scale spatial exploration, spreading and morphogenesis under constraints. Confrontation with experimental data obtained from mycelial thalli of *Podospora anserina*.</footnote>
coordinate to keep track of the current value of time as they do in their definition of $\mathcal{X}$). The flow $\Phi$ describing the deterministic growth process is simply given for any $(\mathfrak{e}, x) \in S$ and $0 \le s \le t$ by

$$ \Phi((\mathfrak{e}, x), s, t) := (\mathfrak{e}, x + \mathfrak{e}(t-s)v). \qquad (2.2) $$

The instantaneous rate at which a division (or reproduction) event happens to an individual with characteristics $(\mathfrak{e}, x)$ is given by

$$ B(\mathfrak{e}, x) := b_1 \mathfrak{e} + b_2 x. \qquad (2.3) $$

The function $B$ is continuous on $S$. Each reproducing individual is replaced by exactly 3 offspring, and so in the notation of [29] we have $p_k(\mathfrak{e}, x) = 1_{\{k=3\}}$ for all $(\mathfrak{e}, x) \in S$. The characteristics at birth of the 3 offspring of an individual of characteristics $z = (\mathfrak{e}, x)$ are given by a triplet $(F_1(z, \Theta), F_2(z, \Theta), F_3(z, \Theta))$, where $\Theta$ is an independent draw from a uniform random variable on $(0, 1)$ and the functions $F_1, F_2, F_3: S \times (0, 1) \to \{0, 1\} \times \mathbb{R}_+$ are given by: for every $(\mathfrak{e}, x) \in S$ and $\theta \in (0, 1)$,

$$ (F_1(\mathfrak{e}, x, \theta), F_2(\mathfrak{e}, x, \theta), F_3(\mathfrak{e}, x, \theta)) \qquad (2.4) $$

$$ := \begin{cases} \left( \left(0, x \frac{\theta(b_1 \mathfrak{e} + b_2 x)}{b_2 x}\right), \left(\mathfrak{e}, x \left(1 - \frac{\theta(b_1 \mathfrak{e} + b_2 x)}{b_2 x}\right)\right), (1, 0) \right) & \text{if } \theta < \frac{b_2 x}{b_1 \mathfrak{e} + b_2 x}, \\ \left((0, x), (1, 0), (1, 0)\right) & \text{if } \theta \ge \frac{b_2 x}{b_1 \mathfrak{e} + b_2 x}. \end{cases} $$

In words, assuming that $\Theta$ follows a uniform distribution on $(0, 1)$ we see that the branching of an open individual ($\mathfrak{e} = 1$) is “lateral” (first line above) with probability $b_2x/(b_1 + b_2x)$, or “apical” with probability $b_1/(b_1 + b_2x)$. Conditionally on $\Theta < b_2x/(b_1 + b_2x)$, the first offspring inherits a fraction $\Theta/[b_2x/(b_1 + b_2x)]$ of the parental length, which indeed corresponds to a uniformly distributed split. If we now consider closed individuals ($\mathfrak{e} = 0$), we have $b_2x/(b_1\mathfrak{e} + b_2x) = 1$ and since we have excluded the case $\theta = 1$, only “lateral” branching can occur and during such an event, the parental length $x$ is split uniformly at random between the first two offspring.

**Remark 2.1.** Note that the function $F_1$ can take the value $(0, 0)$ (which is excluded from $S$), but only when applied to $(1, 0)$. Since open individuals grow at deterministic speed $v > 0$, the amount of time a given individual spends in the state $(1, 0)$ has Lebesgue measure zero. Once we have introduced the Poisson point measure driving the reproduction events below, this will guarantee that the probability that a reproduction event occurs during which an individual with characteristics $(1, 0)$ gives birth to an individual with characteristics $(0, 0)$ is zero. Consequently, provided that the initial state of the population has support in $S$, this property will hold true at any later time with probability one.

On top of the individuals' characteristics, we follow their labels in $\mathcal{U}$ indicating the genealogical relationship between them. Let thus $M_p(\mathcal{U} \times S)$ be the space of all finite point measures on $\mathcal{U} \times S$, equipped with the topology of weak convergence. The state of the population at any time $t$ will take the form

$$ \bar{Z}_t := \sum_{u \in V_t} \delta_{(u, \mathfrak{e}^u, x_t^u)}, \qquad (2.5) $$

where $V_t$ is the index set of all individuals alive at time $t$ (i.e., of all atoms of $\bar{Z}_t$). Restricting our attention to the marginal on $S$ of $\bar{Z}_t$, we shall then obtain the following measure $Z_t \in M_p(S)$:

$$ Z_t := \sum_{u \in V_t} \delta_{(\mathfrak{e}^u, x_t^u)}. \qquad (2.6) $$
The main result of this section is the following theorem. Let $M$ be a Poisson point measure on $\mathbb{R}_+ \times \mathcal{U} \times \mathbb{R}_+ \times (0, 1)$ with intensity $ds \otimes \nu(du) \otimes dz \otimes d\theta$, where $\nu(du)$ denotes the counting measure on $\mathcal{U}$. Let $(\mathcal{F}_t)_{t \ge 0}$ denote the natural filtration associated to $M$. Finally, let $C_b^1(\mathcal{U} \times S)$ stand for the space of all measurable functions on $\mathcal{U} \times S$ that are bounded, continuously differentiable with respect to the variable $x$ and whose first derivatives w.r.t. $x$ are bounded uniformly in $u, \mathfrak{e}$.

**Theorem 2.2.** Let $\bar{Z}^0 \in \mathcal{M}_p(\mathcal{U} \times S)$. Then, there exists a strongly unique $(\mathcal{F}_t)_{t \ge 0}$-adapted càdlàg process $(\bar{Z}_t)_{t \ge 0}$ with values in $\mathcal{M}_p(\mathcal{U} \times S)$ such that $\bar{Z}_0 = \bar{Z}^0$ a.s., and for all $f \in C_b^1(\mathcal{U} \times S)$ and all $t \ge 0$,

$$ \langle \bar{Z}_t, f \rangle = \langle \bar{Z}_0, f \rangle + \int_0^t \int_{\mathcal{U} \times S} v \mathfrak{e} \frac{\partial f}{\partial x}(u, \mathfrak{e}, x) \bar{Z}_s(du, d\mathfrak{e}, dx) ds \quad (2.7) \\ + \int_{[0,t] \times \mathcal{U} \times \mathbb{R}_+ \times (0,1)} \mathbb{1}_{\{u \in V_{s-}, z \le B(\mathfrak{e}^u, x_{s-}^u)\}} \left[ \sum_{i=1}^3 f(\mathfrak{u}i, F_i(\mathfrak{e}^u, x_{s-}^u, \theta)) - f(\mathfrak{u}, \mathfrak{e}^u, x_{s-}^u) \right] M(ds, du, dz, d\theta). $$

We prove Theorem 2.2 and then show that it implies Proposition 1.1.

*Proof of Theorem 2.2.* We only have to check that Assumptions **A** and **B** of Theorem 2.1 in [29] are satisfied.

Concerning **A.1**, we have for all $(\mathfrak{e}, x) \in S$

$$ B(\mathfrak{e}, x) \le (b_1 \vee b_2)(|\mathfrak{e}| + |x|), $$

and so Assumption **A.1** is satisfied with $\gamma = 1$.

Assumption **A.2** is trivially satisfied since for all $(\mathfrak{e}, x) \in S$ and $\theta \in (0, 1)$, we have

$$ \sum_{i=1}^{3} F_i(\mathfrak{e}, x, \theta) \le (2, x) \quad \text{componentwise.} $$

Assumption **A.3** is satisfied since the offspring number is a.s. equal to 3, independently of the parental characteristics.

As concerns Assumption **A.4**, the branching rate of open individuals is bounded from below by $b_1 > 0$, which yields the result in this case. Since elements of $S$ of the form $(0, x)$ satisfy that $x > 0$ by construction, we can write

$$ \int_s^t B(\Phi((0, x), s, r)) dr = b_2 x(t-s) \to +\infty \quad \text{a.s. as } t \to \infty, $$

which is the desired condition.

It remains to check Assumption **B**. Recall from **A.1** that $\gamma = 1$. The infinitesimal generator of the individual trait dynamics corresponding to the flow $\Phi$ is simply given by: for all $f \in C_b^1(S)$ and $(\mathfrak{e}, x) \in S$,

$$ \mathcal{H}f(\mathfrak{e}, x) := v\mathfrak{e} \frac{\partial f}{\partial x}(\mathfrak{e}, x). $$

The function $h: (\mathfrak{e}, x) \mapsto (|\mathfrak{e}| + |x|)^\gamma = \mathfrak{e} + x$ does not belong to the domain of $\mathcal{H}$ since it is not bounded. However, it is easy to construct a sequence $(h_n)_{n \ge 1}$ of functions on $S$ such
that for every $n \ge 1$, $h_n$ and $h$ coincide on the set $\{(\mathfrak{e}, x) \in S : x \le n\}$, $h_n \in C_b^1(S)$ and there exists $c_1, c_2 \ge 0$ such that for all $n \ge 1$ and $(\mathfrak{e}, x) \in S$,

$$ \lim_{n \to \infty} \mathcal{H}h_n(\mathfrak{e}, x) \le c_1(\mathfrak{e} + x) + c_2. $$

For instance, define $h_n(\mathfrak{e}, x) = h(\mathfrak{e}, x)$ if $x \le n$, $h_n(\mathfrak{e}, x) = \mathfrak{e} + (n+1)$ if $x \ge n+1$ and define $h_n(\mathfrak{e}, x)$ for $x \in (n, n+1)$ in such a way that $h_n$ is differentiable with respect to $x$ and its derivative is bounded by 2 on $[n, n+1]$.

Theorem 2.1 in [29], together with Remark 2.1, then yield the result (observing that the martingales $M_{0,\cdot}^f(x)$ in Theorem 2.2 are identically equal to 0 in our case since the individual growth process is deterministic). $\square$

Let us now allow the initial condition $\bar{Z}^0$ to be random (recall that we supposed that $(\Omega, \mathcal{F}, \mathbb{P})$ could be as large as needed), and assume from now on that $\bar{Z}^0$ is such that

$$ \mathbb{E}[\langle \bar{Z}^0, 1 \rangle] < \infty \quad \text{and} \quad \mathbb{E}[\langle \bar{Z}^0, p_l \rangle] < \infty, \qquad (2.8) $$

where we have used again the notation $p_l$ (introduced just before Proposition 1.1) for the projector $(u, \mathfrak{e}, x) \mapsto x$ on the length coordinate.

If we restrict our attention to functions $f$ that are independent of the $\mathcal{U}$-coordinate, Equation (2.7) reads

$$ \langle \mathcal{Z}_t, f \rangle = \langle \mathcal{Z}_0, f \rangle + \int_0^t \int_S v \mathfrak{e} \frac{\partial f}{\partial x}(\mathfrak{e}, x) \mathcal{Z}_s(d\mathfrak{e}, dx) ds \qquad (2.9) \\ + \int_{[0,t] \times \mathcal{U} \times \mathbb{R}_+ \times (0,1)} \mathbb{1}_{\{u \in V_{s-}, z \le B(\mathfrak{e}^u, x_{s-}^u)\}} \left[ \sum_{i=1}^3 f(F_i(\mathfrak{e}^u, x_{s-}^u, \theta)) - f(\mathfrak{e}^u, x_{s-}^u) \right] M(ds, du, dz, d\theta). $$

For $f \equiv 1$, this yields for all $t \ge 0$

$$ \langle \mathcal{Z}_t, 1 \rangle = \langle \mathcal{Z}_0, 1 \rangle + 2 \int_{[0,t] \times \mathcal{U} \times \mathbb{R}_+ \times (0,1)} \mathbb{1}_{\{u \in V_{s-}, z \le B(\mathfrak{e}^u, x_{s-}^u)\}} M(ds, du, dz, d\theta). \qquad (2.10) $$

Taking expectations in the above and using that $B(\mathfrak{e}, x) \le b_1 + b_2x$ for all $(\mathfrak{e}, x) \in S$ and Fubini's theorem, we obtain:

$$ \begin{aligned} \mathbb{E}[\langle \mathcal{Z}_t, 1 \rangle] &= \mathbb{E}[\langle \mathcal{Z}_0, 1 \rangle] + 2 \int_0^t \mathbb{E}[\langle \mathcal{Z}_s, B \rangle] ds \\ &\le \mathbb{E}[\langle \mathcal{Z}_0, 1 \rangle] + 2b_1 \int_0^t \mathbb{E}[\langle \mathcal{Z}_s, 1 \rangle] ds + 2b_2 \int_0^t \mathbb{E}[\langle \mathcal{Z}_s, p_l \rangle] ds \\ &\le \mathbb{E}[\langle \mathcal{Z}_0, 1 \rangle] + 2b_2t \mathbb{E}[\langle \mathcal{Z}_0, p_l \rangle] + (2b_1 + 2b_2vt) \int_0^t \mathbb{E}[\langle \mathcal{Z}_s, 1 \rangle] ds, \end{aligned} \qquad (2.11) $$

where on the last line we have used the fact $(\langle \mathcal{Z}_s, 1 \rangle)_{s \ge 0}$ is a non-decreasing process and therefore

$$ \mathbb{E}[\langle \mathcal{Z}_t, p_l \rangle] \le \mathbb{E}[\langle \mathcal{Z}_0, p_l \rangle] + vt \mathbb{E}[\langle \mathcal{Z}_t, 1 \rangle]. \qquad (2.12) $$

Combining (2.11), Gronwall's lemma and (2.12), we obtain the following bounds.
**Lemma 2.3.** There exists $C_1 > 0$ such that for all $t \ge 0$, we have

$$ \mathbb{E}[\langle \mathcal{Z}_t, 1 \rangle] \le (\mathbb{E}[\langle \mathcal{Z}_0, 1 \rangle] + 2b_2t \mathbb{E}[\langle \mathcal{Z}_0, p_l \rangle]) e^{C_1 t^2} $$

and

$$ \mathbb{E}[\langle \mathcal{Z}_t, p_l \rangle] \le \mathbb{E}[\langle \mathcal{Z}_0, p_l \rangle] + (\mathbb{E}[\langle \mathcal{Z}_0, 1 \rangle] + 2b_2t \mathbb{E}[\langle \mathcal{Z}_0, p_l \rangle]) vte^{C_1 t^2}. $$

Lemma 2.3 gives rather crude bounds on the expectations of the two quantities that control the total branching rate of $(\bar{\mathcal{Z}}_t)_{t \ge 0}$ (or equivalently, of $(\mathcal{Z}_t)_{t \ge 0}$). Its main point is that these expectations are finite at all times. It will be used in the proof of Proposition 1.1.

*Proof of Proposition 1.1.* The proof is rather standard. Recall the notation $F_f$ from (1.4). Conditionning on the value of $\bar{Z}^0$ and using the construction of $(\bar{\mathcal{Z}}_t)_{t \ge 0}$ given in Theorem 2.2, we can write that for every $F \in C_b^1(\mathbb{R})$ and $f \in C_b^1(S)$ (abusing notation and seeing $f$ as a function in $C_b^1(\mathcal{U} \times S)$ independent of the first coordinate), and every $0 \le t \le t'$,

$$ \begin{aligned} & F(\langle \mathcal{Z}_{t'}, f \rangle) - F(\langle \mathcal{Z}_t, f \rangle) - \int_t^{t'} \int_S v \mathfrak{e} F'(\langle \mathcal{Z}_s, f \rangle) \frac{\partial f}{\partial x}(\mathfrak{e}, x) \mathcal{Z}_s(d\mathfrak{e}, dx) ds \\ & = \int_{(t,t'] \times \mathcal{U} \times \mathbb{R}_+ \times (0,1)} \mathbb{1}_{\{u \in V_{s-}, z \le B(\mathfrak{e}^u, x_{s-}^u)\}} \left[ F\left(\langle \mathcal{Z}_{s-}, f \rangle + \sum_{i=1}^3 f(F_i(\mathfrak{e}^u, x_{s-}^u, \theta)) - f(\mathfrak{e}^u, x_{s-}^u)\right) \right. \\ & \qquad \left. - F(\langle \mathcal{Z}_{s-}, f \rangle) \right] M(ds, du, dz, d\theta), \end{aligned} $$

so that

$$ \begin{aligned} & F(\langle \mathcal{Z}_{t'}, f \rangle) - F(\langle \mathcal{Z}_t, f \rangle) - \int_t^{t'} \mathcal{G} F_f(\mathcal{Z}_s) ds \tag{2.13} \\ & = \int_{(t,t'] \times \mathcal{U} \times \mathbb{R}_+ \times (0,1)} \mathbb{1}_{\{u \in V_{s-}, z \le B(\mathfrak{e}^u, x_{s-}^u)\}} \left[ F\left(\langle \mathcal{Z}_{s-}, f \rangle + \sum_{i=1}^3 f(F_i(\mathfrak{e}^u, x_{s-}^u, \theta)) - f(\mathfrak{e}^u, x_{s-}^u)\right) \right. \\ & \qquad \left. - F(\langle \mathcal{Z}_{s-}, f \rangle) \right] M(ds, du, dz, d\theta) \\ & \quad - b_1 \int_t^{t'} \int_S \mathfrak{e} \left\{ F(\langle \mathcal{Z}_{s-}, f \rangle - f(\mathfrak{e}, x) + f(0, x) + 2f(1, 0)) - F(\langle \mathcal{Z}_{s-}, f \rangle) \right\} \mathcal{Z}_{s-}(d\mathfrak{e}, dx) \\ & \quad - b_2 \int_t^{t'} \int_S x \int_0^1 \left\{ F(\langle \mathcal{Z}_{s-}, f \rangle - f(\mathfrak{e}, x) + f(0, (1-\alpha)x) + f(\mathfrak{e}, \alpha x) + f(1, 0)) \right. \\ & \qquad \left. - F(\langle \mathcal{Z}_{s-}, f \rangle) \right\} d\alpha \mathcal{Z}_{s-}(d\mathfrak{e}, dx). \end{aligned} $$

Consequently, using Fubini's theorem we obtain that

$$ \begin{aligned} & \mathbb{E}\left[\left|F(\langle \mathcal{Z}_{t'}, f \rangle) - F(\langle \mathcal{Z}_t, f \rangle) - \int_t^{t'} \mathcal{G} F_f(\mathcal{Z}_s) ds\right|\right] \\ & \le 2 \|F\|_{\infty} \int_t^{t'} \mathbb{E}[\langle \mathcal{Z}_{s-}, B \rangle] ds + 2b_1 \|F\|_{\infty} \int_t^{t'} \mathbb{E}[\langle \mathcal{Z}_{s-}, 1 \rangle] ds \\ & \quad + 2b_2 \|F\|_{\infty} \int_t^{t'} \mathbb{E}[\langle \mathcal{Z}_{s-}, p_l \rangle] ds. \end{aligned} $$
By Lemma 2.3, this quantity is finite and so the quantity on the l.h.s. of (2.13) is integrable for all $0 \le t \le t'$. The martingale property of

$$ \left( F(\langle \mathcal{Z}_t, f \rangle) - F(\langle \mathcal{Z}_0, f \rangle) - \int_0^t \mathcal{G}F_f(\mathcal{Z}_s) ds \right)_{t \ge 0} $$

is then easy to show using the expression given on the r.h.s. of (2.13). Proposition 1.1 is proved. □

### 3 Growth-fragmentation system : Proof of Proposition 1.2

In this section, we analyse the process of mean measures $(n_t)_{t \ge 0}$ defined in (1.9). Namely, we prove Proposition 1.2, which gives conditions under which for every $t \ge 0$, the marginals $n_t^0(dx)$ and $n_t^1(dx)$ of $n_t$ with respect to the type variable are absolutely continuous w.r.t. Lebesgue measure and describes the growth-fragmentation system (1.12) satisfied by the corresponding densities.

Suppose the assumptions of Proposition 1.2 are satisfied. Using the result of Proposition 1.1 with a sequence $(F_j)_{j \ge 1}$ of functions in $C_b^1(\mathbb{R})$ converging to the identity function $\mathrm{Id}$ and whose first derivatives converge to 1 both uniformly over compact subsets of $\mathbb{R}$, together with the dominated convergence theorem (whose use is justified by the bounds obtained in Lemma 2.3), one can write that for every $f \in C_b^1(S)$,

$$ \left( \langle \mathcal{Z}_t, f \rangle - \langle \mathcal{Z}_0, f \rangle - \int_0^t \mathcal{G} \mathrm{Id}_f(\mathcal{Z}_s) ds \right)_{t \ge 0} \qquad (3.1) $$

is a martingale, where we have extended the definition of $\mathcal{G}$ given in (1.5) to write

$$ \mathcal{G} \mathrm{Id}_f(\nu) := \int_S \mathfrak{e}v \frac{\partial f}{\partial x}(\mathfrak{e}, x) \nu(d\mathfrak{e}, dx) + b_1 \int_S \mathfrak{e}\{f(0, x) + 2f(1, 0) - f(\mathfrak{e}, x)\} \nu(d\mathfrak{e}, dx) \\ + b_2 \int_S x \int_0^1 \{f(0, (1-\alpha)x) + f(\mathfrak{e}, \alpha x) + f(1, 0) - f(\mathfrak{e}, x)\} d\alpha \nu(d\mathfrak{e}, dx). $$

Writing that the expectation of the value at any time $t \ge 0$ of the martingale (3.1) is zero, using the decomposition of each $n_t$ introduced in (1.10) and performing a change of variables on the last line, we arrive at

$$ \int_{\mathbb{R}_+} f(1, x)n_t^1(dx) + \int_{\mathbb{R}_+^*} f(0, x)n_t^0(dx) \qquad (3.2) \\ = \int_{\mathbb{R}_+} f(1, x)n_0^1(dx) + \int_{\mathbb{R}_+^*} f(0, x)n_0^1(dx) + v \int_0^t \int_{\mathbb{R}_+} \frac{\partial f}{\partial x}(1, x)n_s^1(dx)ds \\ + 2b_1 f(1, 0) \int_0^t \int_{\mathbb{R}_+} n_s^1(dx)ds + b_1 \int_0^t \int_{\mathbb{R}_+} (f(0, x) - f(1, x))n_s^1(dx)ds \\ + b_2 f(1, 0) \int_0^t \int_{\mathbb{R}_+^*} x(n_s^1(dx) + n_s^0(dx))ds - b_2 \int_0^t \int_{\mathbb{R}_+^*} x(f(1, x)n_s^1(dx) + f(0, x)n_s^0(dx))ds \\ + b_2 \int_0^t \int_{\mathbb{R}_+} \int_0^x (f(1, y) + f(0, y))dy n_s^1(dx)ds + 2b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x f(0, y)dy n_s^0(dx)ds. $$
(Recall from Remark 2.1 that for all $s \ge 0$, $n_s^1(dx)$ gives no mass to $\{0\}$, so that the fifth integral on the r.h.s is well-defined even though $f$ is not defined on $(0,0)$.) This equation will be a key element of the proof below.

Let us first prove that for every $t \ge 0$, the measures $n_t^1$ and $n_t^0$ are absolutely continuous with respect to Lebesgue measure on $\mathbb{R}_+$ and $\mathbb{R}_+^*$, respectively. To do so, we follow the strategy of [34].

Notice first that Eq. (3.2) can be extended to time-dependent test functions $f(t, \mathfrak{e}, x) = f_t(\mathfrak{e}, x)$ defined on $\mathbb{R}_+ \times S$ and of class $C_b^{1,0,1}(\mathbb{R}_+ \times S)$ (i.e., continuously differentiable with respect to the variables $t$ and $x$, bounded and with first derivatives uniformly bounded over $\mathbb{R}_+ \times S$). We have

$$ 
\begin{aligned}
& \int_{\mathbb{R}_+} f_t(1, x)n_t^1(dx) + \int_{\mathbb{R}_+^*} f_t(0, x)n_t^0(dx) && (3.3) \\
&= \int_{\mathbb{R}_+} f_0(1, x)n_0^1(dx) + \int_{\mathbb{R}_+^*} f_0(0, x)n_0^0(dx) + \int_0^t \int_{\mathbb{R}_+} \frac{\partial f_s}{\partial s}(1, x)n_s^1(dx)ds \\
& \quad + \int_0^t \int_{\mathbb{R}_+^*} \frac{\partial f_s}{\partial s}(0, x)n_s^0(dx)ds + v \int_0^t \int_{\mathbb{R}_+} \frac{\partial f_s}{\partial x}(1, x)n_s^1(dx)ds \\
& \quad + 2b_1 \int_0^t f_s(1, 0) \int_{\mathbb{R}_+} n_s^1(dx)ds + b_1 \int_0^t \int_{\mathbb{R}_+^*} (f_s(0, x) - f_s(1, x))n_s^1(dx)ds \\
& \quad + b_2 \int_0^t f_s(1, 0) \int_{\mathbb{R}_+^*} x(n_s^1(dx) + n_s^0(dx))ds \\
& \quad - b_2 \int_0^t \int_{\mathbb{R}_+^*} x(f_s(1, x)n_s^1(dx) + f_s(0, x)n_s^0(dx))ds \\
& \quad + b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x (f_s(1, y) + f_s(0, y))dyn_s^1(dx)ds + 2b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x f_s(0, y)dyn_s^0(dx)ds.
\end{aligned}
$$

We first show that $n_t^1$ is absolutely continuous w.r.t. Lebesgue measure on $\mathbb{R}_+$ for every $t > 0$. We already know from Remark 2.1 that $n_t^1$ has no atom at 0, which will allow us to write all the integrals below over $\mathbb{R}_+^*$ instead of $\mathbb{R}_+$. Let $\phi$ be a nonnegative function in $C_K^1(\mathbb{R}_+)$, the set of all compactly supported functions of class $C^1$ on $\mathbb{R}_+$. Fix $t > 0$. For every $s \in [0, t]$ and $x \in \mathbb{R}_+$, let us define $f_s(x) := \phi(v(t-s) + x)$. It is straightforward to see that $f$ satisfies the following equation:

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial f_s}{\partial s}(x) + v \frac{\partial f_s}{\partial x}(x) = 0, \quad 0 \le s \le t, \\
\\
\displaystyle f_t(x) = \phi(x).
\end{array}
\right. 
\qquad (3.4) 
$$

Let us now set $\varphi_s(\mathfrak{e}, x) = \mathfrak{e}f_s(x)$ for all $s \in [0, t]$ and $(\mathfrak{e}, x) \in S$. Applying (3.3) to $\varphi$, neglecting the negative terms and using (3.4), we can write

$$ 
\begin{aligned}
\int_{\mathbb{R}_+^*} \phi(x)n_t^1(dx) \le & \int_{\mathbb{R}_+^*} \phi(tv + x)n_0^1(dx) + 2b_1 \int_0^t \phi(v(t-s)) \int_{\mathbb{R}_+^*} n_s^1(dx)ds \\
& + b_2 \int_0^t \phi(v(t-s)) \int_{\mathbb{R}_+^*} x(n_s^1(dx) + n_s^0(dx))ds \\
& + b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x \phi(v(t-s) + y)dyn_s^1(dx)ds.
\end{aligned} 
\qquad (3.5) 
$$
Using our assumption that $n_0^1$ has a density with respect to Lebesgue measure and performing the appropriate changes of variables, we can rewrite the inequality (3.5) as

$$ \begin{aligned} \int_{\mathbb{R}_+^*} \phi(x) n_t^1(dx) \le \int_{tv}^\infty \phi(y) n_0^1(y-tv)dy &+ \frac{2b_1}{v} \int_0^{tv} \phi(y) \int_{\mathbb{R}_+^*} n_{t-\frac{y}{v}}^1(dx)dy \\ &+ \frac{b_2}{v} \int_0^{tv} \phi(y) \int_{\mathbb{R}_+^*} x(n_{t-\frac{y}{v}}^1(dx) + n_{t-\frac{y}{v}}^0(dx))dy \\ &+ \frac{b_2}{v} \int_0^{tv} \int_{\mathbb{R}_+^*} \int_0^x \phi(\alpha + y) dy n_{t-\frac{\alpha}{v}}^1(dx)d\alpha. \end{aligned} \quad (3.6) $$

Notice that by the Fubini-Tonneli theorem, we have

$$ \int_0^{tv} \int_{\mathbb{R}_+^*} \int_0^x \phi(\alpha + y) dy n_{t-\frac{\alpha}{v}}^1(dx) d\alpha = \int_0^{tv} \int_{\mathbb{R}_+^*} \phi(\alpha + y) \left( \int_y^\infty n_{t-\frac{\alpha}{v}}^1(dx) \right) dy d\alpha. $$

Finally, plugging the previous expression into (3.6) and performing a last change of variable $(y' = \alpha + y)$, we obtain

$$ \int_{\mathbb{R}_+^*} \phi(x) n_t^1(dx) \le \int_{\mathbb{R}_+^*} \phi(y) H(t, y) dy, \quad (3.7) $$

where

$$ \begin{aligned} H(t, y) := & \mathbb{1}_{\{y>tv\}} n_0^1(y-tv) + \frac{2b_1}{v} \mathbb{1}_{\{y<tv\}} \int_{\mathbb{R}_+^*} n_{t-\frac{y}{v}}^1(dx) \\ & + \frac{b_2}{v} \mathbb{1}_{\{y<tv\}} \int_{\mathbb{R}_+^*} x(n_{t-\frac{y}{v}}^1(dx) + n_{t-\frac{y}{v}}^0(dx)) + \frac{b_2}{v} \int_0^{tv} \int_{y-\alpha}^\infty n_{t-\frac{\alpha}{v}}^1(dx) d\alpha. \end{aligned} $$

Obviously, $H$ is a nonnegative function. Using Lemma 2.3, we can show that $\int_{\mathbb{R}_+^*} H(t, y)dy < \infty$ and since $\phi$ is bounded, the integral on the r.h.s. of (3.7) is also finite. Since this property holds true for every $\phi \in C_K^1(\mathbb{R}_+)$, a simple density argument allows us to conclude that the measure $n_t^1$ is dominated by a measure which is absolutely continuous w.r.t. Lebesgue measure, and this gives us the desired result.

Let us now turn to $n_t^0$ and let $\phi \in C_K^1(\mathbb{R}_+)$. Since individuals of type 0 do not grow, we do not need to consider test functions that depend on time. Instead, we set $f(\mathfrak{e}, x) = (1 - \mathfrak{e})\phi(x)$ for all $(\mathfrak{e}, x) \in S$. Applying (3.2) to $f$ and neglecting the negative terms, we obtain

$$ \begin{aligned} \int_{\mathbb{R}_+^*} \phi(x) n_t^0(dx) \le \int_{\mathbb{R}_+^*} \phi(x) n_0^0(dx) + b_1 \int_0^t \int_{\mathbb{R}_+^*} \phi(x) n_s^1(dx) ds \\ + b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x \phi(y) dy n_s^1(dx) ds + 2b_2 \int_0^t \int_{\mathbb{R}_+^*} \int_0^x \phi(y) dy n_s^0(dx) ds. \end{aligned} $$

Using the fact that $n_0^0$ and all $n_s^1$ have a density with respect to Lebesgue measure, together
with the Fubini-Tonelli theorem, we arrive at

$$ \begin{aligned} \int_{\mathbb{R}_+^*} \phi(x) n_t^0(dx) &\le \int_{\mathbb{R}_+^*} \phi(x) \left( n_0^0(x) + b_1 \int_0^t n_s^1(x) ds + b_2 \int_0^t \int_x^\infty n_s^1(y) dy ds \right. \\ &\qquad \left. + 2b_2 \int_0^t \int_x^\infty n_s^0(dy) ds \right) dx \\ &=: \int_{\mathbb{R}_+^*} \phi(x) G(t, x) dx. \end{aligned} $$

As $G(t, \cdot)$ is nonnegative and integrable, we can use the same arguments as above to conclude that $n_t^0$ is absolutely continuous w.r.t. Lebesgue measure on $\mathbb{R}_+^*$.

The fact that $(n_t^0, n_t^1)_{t \ge 0}$ is a weak solution to the system (1.12) is easily obtained applying Equation (3.2) to test functions such that $f(0, x) = 0$ for all $x > 0$ (to obtain the equation with boundary condition for $n^1$) or such that $f(1, x) = 0$ for all $x \ge 0$ (to obtain the equation for $n^0$). Namely, for $\phi \in C_b^1(\mathbb{R}_+)$ we have

$$ \begin{aligned} \langle n_t^1, \phi \rangle = \langle n_0^1, \phi \rangle &+ v \int_0^t \langle n_s^1, \phi' \rangle ds - \int_0^t \langle n_s^1, (b_1 + b_2(\cdot))\phi \rangle ds + b_2 \int_0^t \langle n_s^1, \int_0^t \phi(y) dy \rangle ds \\ &+ \phi(0) \int_0^t \left( 2b_1 \int_0^\infty n_s^1(y) dy + b_2 \int_0^\infty y(n_s^1(y) + n_s^0(y)) dy \right) ds, \end{aligned} $$

and for $\phi \in C_b^1(\mathbb{R}_+^*)$ we have

$$ \begin{aligned} \langle n_t^0, \phi \rangle = \langle n_0^0, \phi \rangle &+ b_1 \int_0^t \langle n_s^1, \phi \rangle ds - \int_0^t \langle n_s^0, b_2(\cdot)\phi \rangle ds + 2b_2 \int_0^t \langle n_s^0, \int_0^t \phi(y) dy \rangle ds \\ &+ b_2 \int_0^t \langle n_s^1, \int_0^t \phi(y) dy \rangle ds. \end{aligned} $$

This is the weak formulation of (1.12) we aimed for. The proof of Proposition 1.2 is now complete.

# 4 Long time behaviour

In this section we formulate the spectral problems related to the operator in (1.12) and its adjoint (Section 4.1). We then provide the explicit solutions to these problems, the stationary profiles $N_1$ and $N_0$. Finally, in Section 4.2 we prove that the mean measures $n_t^1$ and $n_t^0$ converge towards these profiles as $t \to \infty$ (in a sense that will be made precise in Theorem 4.3). Theorem 1.3 will then be a simple corollary of Theorem 4.3. As explained in Section 1, from now on we suppose that $v = 1$ without loss of generality.

## 4.1 Stationary profiles

For $t > 0$, let us replace $n_t^1$ and $n_t^0$ by $e^{\lambda t} N_1$ and $e^{\lambda t} N_0$, respectively, in (1.12). Doing so, we obtain the following spectral problem related to the operator in (1.12), which captures
the stationary profile of the mean measure:

$$ (N_1)'(x) + (b_1 + b_2x + \lambda)N_1(x) = b_2 \int_x^\infty N_1(y)dy, \quad (4.1a) $$

$$ N_1(0) = 2b_1 \int_0^\infty N_1(y)dy + b_2 \int_0^\infty y(N_1(y) + N_0(y))dy, \quad (4.1b) $$

$$ (b_2x + \lambda)N_0(x) = b_1N_1(x) + 2b_2 \int_x^\infty N_0(y)dy + b_2 \int_x^\infty N_1(y)dy. \quad (4.1c) $$

Somewhat surprisingly, it is possible to find an explicit solution to this system by first solving the equation satisfied by $N_1$ and then solving for $N_0$. See Appendix A.2 for more details.

Next, we shall prove in Appendix A.1 that any eigenvalue $\lambda$ necessarily satisfies

$$ b_1 + \frac{b_2}{\lambda} = \lambda. $$

From now on, we only consider the maximal eigenvalue, which is the unique positive solution to the above equation. It is given by

$$ \lambda = \frac{b_1 + \sqrt{b_1^2 + 4b_2}}{2}. \quad (4.2) $$

Let us now introduce the dual problem. To do so, let us use (1.12) and integration by parts (together with the boundary condition stated in (1.12) to replace $n_t(0)$ by the sum of two integrals) to obtain that for every $f \in C_b^1(S)$, we have

$$ \begin{aligned} \frac{d}{dt} \langle n_t, f \rangle &= \int_{\mathbb{R}_+} \left( (2b_1 + b_2x)f(1,0) + \frac{\partial f}{\partial x}(1,x) - f(1,x)(b_1 + b_2x) + b_2 \int_0^x f(1,y)dy \right. \\ & \qquad \left. + b_1f(0,x) + b_2 \int_0^x f(0,y)dy \right) n_t^1(x)dx \\ & \qquad + \int_{\mathbb{R}_+^*} \left( b_2xf(1,0) - b_2xf(0,x) + 2b_2 \int_0^x f(0,y)dy \right) n_t^0(x)dx. \\ &= \langle n_t, \mathcal{L}f \rangle, \end{aligned} \quad (4.3) $$

where $\mathcal{L}$ is the adjoint operator of the operator acting on $n_t$ in (1.12). The spectral problem associated to $\mathcal{L}$ for the maximal eigenvalue $\lambda$ reads:

$$ \left\{ \begin{array}{l} - \psi_1'(x) + (b_1 + b_2x + \lambda)\psi_1(x) = b_1\psi_0(x) + b_2 \int_0^x \psi_0(y)dy \quad (4.4a) \\ \qquad + b_2 \int_0^x \psi_1(y)dy + \psi_1(0)(2b_1 + b_2x), \end{array} \right. $$

$$ (b_2x + \lambda)\psi_0(x) = 2b_2 \int_0^x \psi_0(y)dy + b_2x\psi_1(0), \quad (4.4b) $$

The eigenvector $(\psi_1, \psi_0)$ will allow us to quantify the influence of the initial condition on the growth of the population size.

Here again, the spectral problem can be solved to obtain an explicit expression for $\psi_0$ and $\psi_1$. This leads to the following result.
**Proposition 4.1.** The following quadruplet of non-negative functions $(N_0, N_1, \psi_0, \psi_1)$ is solution to the spectral problem given by (4.1) and (4.4) :

$$
\begin{align*}
N_1(x) &= (b_1 + b_2x + \lambda)e^{-\int_0^x (b_1+b_2y+\lambda)dy}, \\
N_0(x) &= \frac{N_1(x)}{(b_2x + \lambda)^2}(b_2 + b_1(b_2x + \lambda)) + \frac{b_2e^{-\int_0^x (b_1+b_2y+\lambda)dy}}{(b_2x + \lambda)^3}(2b_2 + b_1(b_2x + \lambda)), \\
\psi_1(x) &= c_0\left(1 + \frac{b_2}{\lambda}x\right), \\
\psi_0(x) &= c_0\frac{b_2x}{\lambda},
\end{align*}
$$

where $c_0 = \frac{\lambda^2}{\lambda^2+b_2}$. Besides, this solution satisfies

$$
\int_0^\infty (N_1(x) + N_0(x)) dx = 2, \quad \text{and} \quad \int_0^\infty (\psi_1(y)N_1(y) + \psi_0(y)N_0(y)) dy = 1. \quad (4.5)
$$

Note that the definition of $\psi_1$, $\psi_0$ in Proposition 4.1 is equivalent to the definition of $\psi$
given in (1.16) in that $\psi(\mathfrak{e}, x) = \psi_{\mathfrak{e}}(x)$.

The proof of Proposition 4.1 is given in Appendix A.2. The quadruplet of functions
provides the unique solution to the spectral problem associated to the maximal eigenvalue
$\lambda$ under the normalising conditions (4.5). Uniqueness in the weighted space related to a
particular Lyapunov function $V_\gamma$ will be obtained in the next section, see Theorem 4.3.

**4.2 Convergence of the mean measure**

In this section, we prove Theorem 1.3 about the long-time convergence of the mean density
of open and closed individuals (or, in view of our application, of external and internal
filaments). In fact, we shall prove a more general statement involving the semigroup
associated to the process ($\mathcal{Z}_t$)$_{t \ge 0}$, of which Theorem 1.3 will be an easy corollary. Before
we present this result, we need some notational preparation. The approach we shall adopt
is the semigroup approach of [7].

Recall the definition of the function $\psi$ given in (1.16), motivated by the result of Proposition 4.1. For any $\gamma \ge 2$ and for all $(\mathfrak{e}, x) \in S$, let us define

$$
V_{\gamma}(\mathfrak{e}, x) = \psi(\mathfrak{e}, x) + x^{\gamma} + 1. \tag{4.6}
$$

Observe that when $\gamma = 2$, we recover the definition of $V$ given in (1.16). Let $\mathcal{B}(V_\gamma)$ denote the set of all measurable functions $f : S \to \mathbb{R}$ such that the following quantity is finite:

$$
\|f\|_{\mathcal{B}(V_\gamma)} := \sup_{z \in S} \frac{|f(z)|}{V_\gamma(z)}. \qquad (4.7)
$$

Let $\mathcal{M}(V_\gamma)$ denote the set of all signed measures on $S$ that integrate $V_\gamma$. The space $\mathcal{M}(V_\gamma)$ is endowed with the weighted total variation norm

$$
\|\mu\|_{\mathcal{M}(V_\gamma)} := \sup_{\|f\|_{\mathcal{B}(V_\gamma)} \le 1} \left| \int_S f(z) \mu(dz) \right|. \quad (4.8)
$$

By Proposition 1.1, the stochastic process ($\mathcal{Z}_t$)$_{t \ge 0}$ is well-defined for any initial condition made of a single atom at some $z = (\mathfrak{e}, x) \in S$. We can thus define, for any $t \ge 0$ and any nonnegative measurable function $f$ on $S$ :

$$
M_t f(z) = \mathbb{E}_{\delta_z}[\langle \mathcal{Z}_t, f \rangle] \in [0, +\infty]. \quad (4.9)
$$
Let us set

$$ \mathfrak{B} = \bigcup_{\gamma \ge 2} \mathcal{B}(V_{\gamma}). $$

The following result extends $M = (M_t)_{t \ge 0}$ to this set of functions, on which it takes values in the set of finite functions on $S$ and satisfies the semigroup property.

**Lemma 4.2.** (i) For any $\gamma \ge 2$, there exists $C_\gamma > 0$ such that for all $t \ge 0$,

$$ M_t V_\gamma(z) \le e^{C_\gamma t} V_\gamma(z), \quad \forall z \in S. $$

(ii) For any nonnegative $f \in \mathfrak{B}$, $M_t f(z)$ is finite for all $z \in S$. We can therefore extend the definition of $M$ to $\mathfrak{B}$ as follows: For any $\gamma \ge 2$, $f \in \mathcal{B}(V_\gamma)$ and $t \ge 0$, we set

$$ M_t f(z) = \mathbb{E}_{\delta_z} [\langle \mathcal{Z}_t, f \rangle] := M_t f_+(z) - M_t f_-(z), \quad \forall z \in S, $$

where $f_+$ (resp. $f_-$) is the positive (resp. negative part) of $f$. We have $M_t f \in \mathcal{B}(V_\gamma)$.

(iii) $(M_t)_{t \ge 0}$ is a positive semigroup on $\mathfrak{B}$ and satisfies

$$ M_t f(\mathfrak{e}, x) = f(\mathfrak{e}, x + \mathfrak{e}t)e^{-\int_0^t (b_1 \mathfrak{e} + b_2(x+\mathfrak{e}s)) ds} \\ + \int_0^t e^{-\int_0^s (b_1 \mathfrak{e} + b_2(x+\mathfrak{e}s')) ds'} \int_S M_{t-s} f(\bar{\mathfrak{e}}, \bar{x}) Q(\mathfrak{e}, x + \mathfrak{e}s, d\bar{\mathfrak{e}}, d\bar{x}) ds, $$

where for $(\mathfrak{e}, x) \in S$, we have

$$ Q(\mathfrak{e}, x, d\bar{\mathfrak{e}}, d\bar{x}) = \mathbb{1}_{\{\bar{x} \le x\}} \left[ b_1 \mathfrak{e} \left[ 2\delta_1(d\bar{\mathfrak{e}}) \otimes \delta_0(d\bar{x}) + \delta_0(d\bar{\mathfrak{e}}) \otimes \delta_x(d\bar{x}) \right] \right. \\ \left. + b_2 x \left[ \delta_1(d\bar{\mathfrak{e}}) \otimes \delta_0(d\bar{x}) + \delta_0(d\bar{\mathfrak{e}}) \otimes \frac{d\bar{x}}{x} + \delta_\mathfrak{e}(d\bar{\mathfrak{e}}) \otimes \frac{d\bar{x}}{x} \right] \right]. $$

*Proof. (Proof of Lemma 4.2)* To prove (i), let us observe that the adjoint operator $\mathcal{L}$ introduced in (4.3) can be rewritten in a more compact way as follows:

$$ \mathcal{L} f(\mathfrak{e}, x) := \mathfrak{e} \frac{\partial f}{\partial x} + b_1 \mathfrak{e} (f(0, x) + 2f(1, 0) - f(\mathfrak{e}, x)) \\ + b_2 x \int_0^1 (f(0, (1-\alpha)x) + f(\mathfrak{e}, \alpha x) + f(1, 0) - f(\mathfrak{e}, x)) d\alpha, \quad (4.10) $$

$$ = \mathfrak{e} \frac{\partial f}{\partial x} - (b_1 \mathfrak{e} + b_2 x) f(\mathfrak{e}, x) + \int_S f(\bar{\mathfrak{e}}, \bar{x}) Q(\mathfrak{e}, x, d\bar{\mathfrak{e}}, d\bar{x}), \quad (4.11) $$

and that if we write (1.5) with $F = \text{Id}$, we have for all $\nu \in \mathcal{M}_p(S)$

$$ \mathcal{G} \text{Id}_f(\nu) = \int_S \mathcal{L} f(\mathfrak{e}, x) \nu(d\mathfrak{e}, dx). \quad (4.12) $$

For convenience, let us define the function $h_\gamma : (\mathfrak{e}, x) \mapsto x^\gamma + 1$, so that $V_\gamma = \psi + h_\gamma$. Since $\psi$ is an eigenfunction for $\mathcal{L}$ associated to the eigenvalue $\lambda$, we have $\mathcal{L}\psi = \lambda\psi$. Furthermore, using (4.11) we can write that for every $(\mathfrak{e}, x) \in S$,

$$ \begin{aligned} \mathcal{L} h_\gamma(\mathfrak{e}, x) &= \mathfrak{e} \gamma x^{\gamma-1} - (b_1 \mathfrak{e} + b_2 x)(x^\gamma + 1) + 2b_1 \mathfrak{e} \times 1 + b_1 \mathfrak{e}(x^\gamma + 1) + b_2 x \times 1 \\ &\quad + 2b_2 \int_0^x (\bar{x}^\gamma + 1) d\bar{x} \\ &= \mathfrak{e} \gamma x^{\gamma-1} - b_2 x^{\gamma+1} + 2b_1 \mathfrak{e} + \frac{2b_2}{\gamma+1} x^{\gamma+1} + 2b_2 x \\ &= \mathfrak{e} \gamma x^{\gamma-1} - b_2 \frac{\gamma-1}{\gamma+1} x^{\gamma+1} + 2b_1 \mathfrak{e} + 2b_2 x. \end{aligned} \quad (4.13) $$
Combining these two results, we obtain that for every $(\mathfrak{e}, x) \in S$,

$$ \mathcal{L}V_{\gamma}(\mathfrak{e}, x) = \lambda\psi(\mathfrak{e}, x) + \mathfrak{e}\gamma x^{\gamma-1} - b_2 \frac{\gamma-1}{\gamma+1} x^{\gamma+1} + 2b_1 \mathfrak{e} + 2b_2 x. \quad (4.14) $$

Now recall from Proposition 4.1 (or, equivalently, from Equation (1.16)) that

$$ \psi(\mathfrak{e}, x) = c_0 \left( \mathfrak{e} + \frac{b_2}{\lambda} x \right) \quad \text{with } c_0 = \frac{\lambda^2}{\lambda^2 + b_2}, $$

so that

$$ 2b_2x \le 2\lambda \left( \mathfrak{e} + \frac{b_2}{\lambda} x \right) \le \frac{2\lambda}{c_0} \psi(\mathfrak{e}, x) \quad \text{and} \quad \lambda\psi \le \frac{\lambda}{c_0} \psi. \quad (4.15) $$

Likewise, we have

$$ 2b_1 \mathfrak{e} \le \frac{2b_1}{c_0} \psi(\mathfrak{e}, x). \quad (4.16) $$

Next, since $\mathfrak{e} \le 1$, we can write that

$$ \mathfrak{e}\gamma x^{\gamma-1} - b_2 \frac{\gamma-1}{\gamma+1} x^{\gamma+1} \le \mathfrak{e} \left( \gamma x^{\gamma-1} - b_2 \frac{\gamma-1}{\gamma+1} x^{\gamma+1} \right). $$

When $\mathfrak{e} = 0$ the bound on the r.h.s. is zero, while when $\mathfrak{e} = 1$ the expression on the r.h.s. is bounded by some constant $x_0 = x_0(\gamma) > 0$. Consequently, the quantity on the l.h.s. is bounded by $\mathfrak{e}x_0 \le (x_0/c_0)\psi$. Combining the above, we obtain that

$$ \mathcal{L}V_{\gamma}(\mathfrak{e}, x) \le \frac{3\lambda + 2b_1 + x_0}{c_0} \psi(\mathfrak{e}, x). \quad (4.17) $$

Using the Kolmogorov equation (valid for $f \in C_b^1(S)$)

$$ \frac{d}{dt} \mathbb{E}[\langle \mathcal{Z}_t, f \rangle] = \mathbb{E}[\langle \mathcal{Z}_t, \mathcal{L}f \rangle]. $$

with a sequence $(f_n)_{n \ge 0}$ of functions increasing to $V_{\gamma}$, together with a standard monotone convergence argument, the fact that $\psi \le V_{\gamma}$ and finally Gronwall's lemma, we obtain that for any $z \in S$ and $t \ge 0$,

$$ M_t V_{\gamma}(z) \le e^{(3\lambda+2b_1+x_0)t/c_0} V_{\gamma}(z). \quad (4.18) $$

(ii) is a direct consequence of (i) noticing that the definition of $M$ on the embedded sets $\mathcal{B}(V_{\gamma})$ is compatible since it coincides with $\mathbb{E}[\langle \mathcal{Z}_t, f \rangle | \mathcal{Z}_0 = \delta_z]$. The semigroup property is classical : it is a consequence of the branching Markov property of $\mathcal{Z}$. Finally, the proof of Duhamel's formula in (iii) comes as usual by conditioning on the first jump of $\mathcal{Z}$ and using the strong Markov property. $\square$

We are ready to state the convergence result in its full generality. For every bounded measurable $f : S \to \mathbb{R}$, we define the following measure $N$ on $S$, in the same spirit as the decomposition (1.10):

$$ \langle N, f \rangle = \int_{\mathbb{R}_+} f(1, x) N_1(dx) + \int_{\mathbb{R}_+^*} f(0, x) N_0(dx). \quad (4.19) $$

The fact that we know (explicitly here) an eigenfunction $\psi$ allows us to invoke a Doob $h$-transform. This method is powerful to study non-conservative semigroups [31], and in
particular to study the first moment properties of branching processes. We can then derive ergodic estimates for our semigroup $M$ from a Harris ergodic theorem applied to the associated conservative semigroup, as recently achieved in [14] for other growth-fragmentation PDEs. Instead, here we obtain these results directly by applying Theorem 2.1 in [7]. The two methods are equivalent but the latter is more convenient in our framework. It also allows the extension of the results to models where the positive eigenfunction is not known a priori, which would be the case when the growth rate or fragmentation are different.

**Theorem 4.3.** Let $\gamma \ge 2$. There exist $C, w > 0$, depending on $\gamma$, such that for all $t \ge 0$ and $\mu \in \mathcal{M}(V_\gamma)$,

$$ \|e^{-\lambda t}\mu M_t - \langle \mu, \psi \rangle N\|_{\mathcal{M}(V_\gamma)} \le C e^{-wt}\|\mu\|_{\mathcal{M}(V_\gamma)}. $$

Let us mention that combining the proof below and the results of [7] provides a lower bound on the spectral gap and speed of convergence $w$ in terms of the parameters of the growth-fragmentation model, see forthcoming Remark 4.4.

Before we prove Theorem 4.3, let us show how we can deduce Theorem 1.3 from it.

*Proof. (Proof of Theorem 1.3)* Take $\mu = n_0$ (where $n_0$ is the finite measure on $S$ constructed from the densities $n_0^0$ and $n_0^1$ as in (1.10)) and observe that the function $\varphi \equiv 1$ belongs to $\mathcal{B}(V_\gamma)$ for any $\gamma \ge 2$ and satisfies $\|\varphi\|_{\mathcal{B}(V_\gamma)} \le 1$. Since

$$ \|n_0\|_{\mathcal{M}(V_\gamma)} = \sup_{\|f\|_{\mathcal{B}(V_\gamma)} \le 1} \left| \int_S f(z)n_0(dz) \right| \le \int_S V_\gamma(z)n_0(dz), $$

taking $\gamma = 2$ we obtain that (1.17) is indeed satisfied. $\square$

*Proof. (Proof of Theorem 4.3)* We write $f \lesssim g$ to mean that there exists a constant $C > 0$ (also independent of the time variable when this notation is used with $M_t f$ over a finite time interval $[0, T]$) such that $f \le Cg$. For some $R_\gamma > 0$ that will be chosen later, let us define

$$ \mathbb{K}_\gamma = \{z \in S : V_\gamma(z) \le R_\gamma \psi(z)\} $$

and first observe that $\psi \le V_\gamma$ on $S$ and $V_\gamma \lesssim \psi$ on $\mathbb{K}_\gamma$. Second, using (4.18) and the fact that $\psi$ is an eigenfunction, we obtain that for every $T > 0$,

$$ MV_\gamma \lesssim V_\gamma \quad \text{and} \quad M\psi \gtrsim \psi \quad \text{on} \quad [0, T] \times S. $$

To apply Theorem 2.1 in [7] and obtain the desired estimate, we need to check the following assumption.

**Assumption A.** There exist $\tau > 0$, $\beta > \alpha > 0$, $\theta \ge 0$, $(c, d) \in (0, 1]^2$, $\mathbb{K} \subset S$ and $\nu$ a probability measure on $S$ supported in $\mathbb{K}$ such that

$$ (A1) \quad M_\tau V_\gamma \le \alpha V_\gamma + \theta 1_{\mathbb{K}} \psi, $$

$$ (A2) \quad M_\tau \psi \ge \beta \psi, $$

$$ (A3) \quad \text{For all } z \in \mathbb{K} \text{ and all nonnegative function } f \in \mathcal{B}(V_\gamma/\psi), $$

$$ M_\tau(f\psi)(z) \ge c \langle \nu, f \rangle M_\tau \psi(z), $$
(A4) For every integer $n \ge 1$,

$$d \sup_{z \in \mathbb{K}} \frac{M_{n\tau}\psi(z)}{\psi(z)} \le \langle \nu, \frac{M_{n\tau}\psi}{\psi} \rangle.$$

Recall that in view of (3.1) and (4.12), we have for every $f \in \mathcal{B}(V_\gamma)$

$$\frac{\partial(M_t f)}{\partial t} = M_t(\mathcal{L}f). \qquad (4.20)$$

Informally, Assumption (A1) corresponds to a Lyapunov-type condition to control the trait distribution in the population (and to show that these traits tend to be confined in compact sets). Assumption (A3) corresponds to a Doeblin (small set)-type condition and guarantees a local mixing property which, in turn, yields the exponential convergence of the trait distribution. The additional Assumptions (A2) and (A4) are needed to exploit these estimates and to control the original (non-conservative) semigroup $M$. The fact that we use the eigenfunction $\psi$ there makes them particularly natural to check.

Let us first observe that for any $\tau > 0$, (A2) and (A4) are satisfied with $\beta = e^{\lambda\tau}$ and $d = 1$. Indeed, $\psi$ is an eigenfunction of the adjoint operator, so that $\mathcal{L}\psi = \lambda\psi$. Consequently, we have for any $t \ge 0$,

$$M_t \psi = e^{\lambda t} \psi.$$

The rest of the proof consists of two steps. In **Step 1** we check (A1), while **Step 2** is devoted to checking (A3). Suitable values for $R_\gamma$ will be fixed in **Step 1**, while the value of $\tau$ will be fixed in **Step 2**.

**Step 1.** Using (4.17), we can write that for every $(\mathfrak{e}, x) \in S$,

$$\begin{aligned} M_t V_\gamma(\mathfrak{e}, x) &\le V_\gamma(\mathfrak{e}, x) + \frac{3\lambda + 2b_1 + x_0}{c_0} \int_0^t M_s \psi(\mathfrak{e}, x) ds \\ &= V_\gamma(\mathfrak{e}, x) + \frac{(3\lambda + 2b_1 + x_0)(e^{\lambda t} - 1)}{c_0} \frac{1}{\lambda} \psi(\mathfrak{e}, x). \end{aligned}$$

Now, using that on $\mathbb{K}_\gamma^c$ we have $\psi \le \frac{1}{R_\gamma} V_\gamma$, we can write

$$\begin{aligned} M_t V_\gamma(\mathfrak{e}, x) &\le V_\gamma(\mathfrak{e}, x) + \frac{(3\lambda + 2b_1 + x_0)(e^{\lambda t} - 1)}{\lambda c_0} \frac{1}{R_\gamma} V_\gamma(\mathfrak{e}, x) \\ &\quad + \frac{(3\lambda + 2b_1 + x_0)(e^{\lambda t} - 1)}{\lambda c_0} \psi(\mathfrak{e}, x) \mathbb{1}_{\mathbb{K}_\gamma}(\mathfrak{e}, x). \end{aligned}$$

Now, suppose we have fixed a value for $\tau$ (which we shall do in the next step). To obtain (A1), it is natural to set

$$\alpha := 1 + \frac{(3\lambda + 2b_1 + x_0)(e^{\lambda\tau} - 1)}{\lambda c_0} \frac{1}{R_\gamma} \qquad (4.21)$$

and

$$\theta := \frac{(3\lambda + 2b_1 + x_0)(e^{\lambda\tau} - 1)}{\lambda c_0}. \qquad (4.22)$$
To obtain that $\alpha < \beta = e^{\lambda\tau}$, it suffices to take $R_{\gamma}$ large enough. Note that an appropriate lower bound on $R_{\gamma}$ ensuring that the latter condition is satisfied depends on the constant multiplying $e^{\lambda\tau}$ in (4.21), but can be taken to be independent of $\tau$.

Finally, with our choice of $R_{\gamma}$ we should justify that $\mathbb{K}_{\gamma}$ is non empty and bounded. As we can choose $R_{\gamma} > C_0 + 2$, it is easy to see that $\mathbb{K}_{\gamma}$ then contains all points $(1, x)$ such that $0 \le x \le 1$. Moreover, since for any $\mathfrak{e} \in \{0, 1\}$ we have

$$ \lim_{x \to \infty} \frac{\psi(\mathfrak{e}, x)}{x^{\gamma}} = 0, $$

the values of $x$ for which $(\mathfrak{e}, x)$ belongs to $\mathbb{K}_{\gamma}$ must be bounded. As a consequence, there exist $k_0, K_0$ and $K_1$ such that

$$ \mathbb{K}_{\gamma} = (\{0\} \times [k_0, K_0]) \cup (\{1\} \times [0, K_1]). $$

Pairs of the form $(0, x)$ where $x$ is small enough do not belong to $\mathbb{K}_{\gamma}$ as $V_{\gamma}(0, x) \ge 1$ while $\psi(0, x)$ tends to 0 as $x$ tends to 0, which means that $k_0 > 0$ in the above.

**Step 2.** Thanks to the properties of $\psi$, we can reformulate Assumption (A3) in the following way: there exist $\tau > 0$, $c \in (0, 1]$, and $\nu$ a probability measure on $S$ supported in $\mathbb{K}_{\gamma}$ such that for all $z \in \mathbb{K}_{\gamma}$ and all nonnegative function $f \in \mathcal{B}(V_{\gamma}/\psi)$, we have

$$ \frac{M_{\tau}(f\psi)(z)}{e^{\lambda\tau}\psi(z)} \ge c\langle\nu, f\rangle. $$

To start with some intuitive argument, observe that the measure $\nu$ is meant to put its weight on a compact subset of $S$ to which the lines of descent of “typical” individuals sampled from the population come back recurrently. Because of the growth-fragmentation dynamics, in which “lateral” branching/fragmentation splits an individuals’ length uniformly at random, we expect an appropriate measure $\nu$ to be absolutely continuous with respect to Lebesgue measure on the “length” coordinate and to have support in a connected set of the form $(\{0\} \times [a_0, b_0]) \cup (\{1\} \times [0, b_1])$ with $0 < a_0 < b_0$ and $b_1 > 0$ (recall that $(0, 0) \notin S$). This is indeed what we shall obtain in (4.36).

To rigorously construct these objects, let us define a new semigroup $P = (P_t)_{t \ge 0}$ as follows. For every $z \in S$, $t \ge 0$ and every nonnegative $f \in \mathcal{B}(V_{\gamma}/\psi)$, let us set

$$ P_t f(z) := \frac{M_t(\psi f)(z)}{e^{\lambda t} \psi(z)}. \qquad (4.23) $$

Defined in this way, $P$ is a conservative semigroup (take $f = 1$ to see the conservation of mass property). The above condition can thus be rewritten: there exist $\tau > 0$, $c \in (0, 1]$, and $\nu$ a probability measure on $S$ supported in $\mathbb{K}_{\gamma}$ such that for all $z \in \mathbb{K}_{\gamma}$ and all nonnegative $f \in \mathcal{B}(V_{\gamma}/\psi)$, we have

$$ P_{\tau} f(z) \ge c \int_{\mathbb{K}_{\gamma}} f(z') \nu(dz'). \qquad (4.24) $$

In order to prove the above statement, we shall first analyse the generator $\mathcal{A}$ of $P$. For an appropriate test function $f$, we compute

$$ \mathcal{A}f(\mathfrak{e}, x) = \frac{\partial(P_t f)}{\partial t}\Big|_{t=0} (\mathfrak{e}, x). $$
Using the formulation of $\mathcal{L}$ given (4.11) to pass from the first to the second line, and
then the fact that $\lambda = (\mathcal{L}\psi)/\psi$ to pass from the second to the third line, we can write

$$
\begin{align*}
\left. \frac{\partial}{\partial t} M_t(\psi f) \right|_{t=0} (\mathfrak{e}, x) &= \frac{\mathcal{L}(\psi f)}{\psi}(\mathfrak{e}, x) - \lambda f(\mathfrak{e}, x) \\
&= \mathfrak{e} \frac{\partial \psi}{\partial x}(\mathfrak{e}, x) \frac{f(\mathfrak{e}, x)}{\psi(\mathfrak{e}, x)} + \mathfrak{e} \frac{\partial f}{\partial x}(\mathfrak{e}, x) - (b_1 \mathfrak{e} + b_2 x) f(\mathfrak{e}, x) \\
&\quad + \int_S \frac{f(\mathfrak{e}', x') \psi(\mathfrak{e}', x')}{\psi(\mathfrak{e}, x)} Q(\mathfrak{e}, x, d\mathfrak{e}', dx') - \lambda f(\mathfrak{e}, x) \\
&= \mathfrak{e} \frac{\partial f}{\partial x}(\mathfrak{e}, x) + \int_S (f(\mathfrak{e}', x') - f(\mathfrak{e}, x)) \frac{\psi(\mathfrak{e}', x')}{\psi(\mathfrak{e}, x)} Q(\mathfrak{e}, x, d\mathfrak{e}', dx').
\end{align*}
$$

Let us define

$$
\bar{B}(\mathfrak{e}, x) := \int_S \frac{\psi(\mathfrak{e}', x')}{\psi(\mathfrak{e}, x)} Q(\mathfrak{e}, x, d\mathfrak{e}', dx'), \quad \forall (\mathfrak{e}, x) \in S. \quad (4.25)
$$

Then, as before, the semigroup $(P_t)_{t \ge 0}$ admits the following representation:

$$
P_t f(\mathfrak{e}, x) = f(\mathfrak{e}, x + \mathfrak{e}t) e^{-\int_0^t \bar{B}(\mathfrak{e}, x + \mathfrak{e}s) ds} \quad (4.26)
$$

$$
+ \int_0^t e^{-\int_0^s \bar{B}(\mathfrak{e}, x + \mathfrak{e}s') ds'} \int_S P_{t-s} f(\mathfrak{e}', x') \bar{Q}(\mathfrak{e}, x + \mathfrak{e}s, d\mathfrak{e}', dx') ds, \quad (4.27)
$$

where for $(\mathfrak{e}, x) \in S$, we have set

$$
\bar{Q}(\mathfrak{e}, x, d\mathfrak{e}', dx') = \frac{\psi(\mathfrak{e}', x')}{\psi(\mathfrak{e}, x)} Q(\mathfrak{e}, x, d\mathfrak{e}', dx').
$$

Here, $\psi$ should be understood as a weight function.

Before we start building the measure $\nu$, we first compute the death rate $\bar{B}$. Using
(4.25), we readily obtain that for $(\mathfrak{e}, x) \in S$,

$$
\bar{B}(\mathfrak{e}, x) = \frac{1}{\psi(\mathfrak{e}, x)} \left( 2b_1 \mathfrak{e} + b_1 \mathfrak{e} \frac{b_2 x}{\lambda} + b_2 x + \frac{b_2^2 x^2}{\lambda} + b_2 x \mathfrak{e} \right).
$$

Let us first remark that $\bar{B}(\mathfrak{e}, x) > 0$ for all $(\mathfrak{e}, x) \in S$. In addition, it is straightforward to check that there exist $\varrho_1, \varrho_2 > 0$ depending only on $b_1, b_2$ and $\lambda$ such that for every $A > 0$, we have

$$
\sup_{\mathfrak{e} \in \{0,1\}, x \le A} \bar{B}(\mathfrak{e}, x) \le \varrho_1 + \varrho_2 A. \quad (4.28)
$$

The latter will be very useful when proving (4.24), as for $x \le K_0 \vee K_1$ we shall have

$$
e^{-\int_{0}^{t} \bar{B}(\mathfrak{e}, x + \mathfrak{e}s') ds'} \geq e^{-t(\varrho_{1} + \varrho_{2}(K_{0} \vee K_{1} + t))}. \quad (4.29)
$$

Let us now start from an individual $(1, x) \in \mathbb{K}_\gamma$ and fix $t > 2(K_0 \vee K_1)$. From (4.27), keeping only the two open filaments of size zero coming from the term describing the apical branching, we have

$$
P_t f(1, x) \ge \int_0^t e^{-\int_0^s \bar{B}(1, x+s') ds'} 2(P_{t-s} f)(1, 0) \frac{1}{1 + \frac{b_2}{\lambda}(x+s)} ds.
$$
Now, for one of the terms of the form $(P_{t-s}f)(1,0)$, we only keep the growth term (4.26). For the other one, we only keep the term corresponding to apical branching and an individual of type 0 coming from it. Doing so, we obtain

$$P_t f(1, x) \geq H_t + J_t, \quad (4.30)$$

where

$$H_t = b_1 \int_0^t e^{-\int_0^s \bar{B}(1, x+s')ds'} e^{-\int_0^{t-s} \bar{B}(1, s')ds'} f(1, s) \frac{1}{1 + \frac{b_2}{\lambda}(x+s)} ds, \\ J_t = b_1^2 \int_0^t e^{-\int_0^s \bar{B}(1, x+s')ds'} \frac{1}{1 + \frac{b_2}{\lambda}(x+s)} \int_0^{t-s} e^{-\int_0^u \bar{B}(1, u')du'} P_{t-s-u}f(0, u) \frac{\psi(0, u)}{\psi(1, u)} du ds.$$

Now, for the first term we use that $x \leq K_1$ and $t \geq K_1$. It comes, using (4.29)

$$H_t \geq \frac{b_1 e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{1 + \frac{b_2}{\lambda}(K_1 + t)} \int_0^{K_1} f(1, s) ds.$$

For the second term, we find yet another lower bound, by keeping only the growth term (4.26) and apply Fubini's theorem. It comes, after using that $x \leq K_1$ and $t \geq 2K_0$,

$$J_t \geq \frac{b_1^2}{1 + \frac{b_2}{\lambda}(K_1 + t)} \int_0^t e^{-\int_0^s \bar{B}(1, x+s')ds'} \int_0^{t-s} e^{-\int_0^u \bar{B}(1, u')du'} e^{-(t-s-u)\bar{B}(0, u)} f(0, u) \times \frac{(b_2 u)/\lambda}{1 + (b_2 u)/\lambda} du ds.$$

Using again (4.29) and defining

$$C_1(t) := \frac{b_1^2 \frac{b_2}{\lambda} e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{(1 + \frac{b_2}{\lambda}(K_1 + t))(1 + \frac{b_2}{\lambda}t)},$$

we obtain

$$J_t \geq C_1(t) \int_0^t \int_0^{t-s} uf(0, u)du ds \geq C_1(t) \int_0^t uf(0, u)(t-u)du \\ \geq K_0 C_1(t) \int_0^{K_0} uf(0, u)du.$$

Plugging the bounds for $H_t$ and $J_t$ into (4.30), we obtain

$$P_t f(1, x) \geq C_2(t) \int_0^{K_1} f(1, s) ds + K_0 C_1(t) \int_0^{K_0} uf(0, u) du, \quad (4.31)$$

where

$$C_2(t) := \frac{b_1 e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{1 + \frac{b_2}{\lambda}(K_1 + t)}.$$

Hence,

$$P_t f(1, x) \geq C_2(t) \int_0^{K_1} f(1, s) ds + K_0 C_1(t) k_0 \int_{k_0}^{K_0} f(0, u) du. \quad (4.32)$$
Let us now start from $(0, x) \in \mathbb{K}_{\gamma}$. Keeping only the term corresponding to lateral branching in (4.27) yields

$$P_t f(0, x) \geq \int_0^t e^{-\bar{B}(0,x)s} P_{t-s} f(1, 0) \frac{\psi(1, 0)}{\psi(0, x)} b_2 x ds = \lambda \int_0^t e^{-\bar{B}(0,x)s} P_{t-s} f(1, 0) ds.$$

Now, using (4.31) with $t$ replaced by $t-s$ and $x=0$, we obtain

$$P_t f(0, x) \geq \lambda \int_0^t e^{-\bar{B}(0,x)s} \left( C_2(t-s) \int_0^{K_1} f(1, u) du + K_0 C_1(t-s) \int_0^{K_0} u f(0, u) du \right) ds.$$

Using again (4.29) together with the definition of the functions $C_1$ and $C_2$, we obtain

$$P_t f(0, x) \geq \lambda \frac{b_1 e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{1 + \frac{b_2}{\lambda}(K_1 + t)} \int_0^t \int_0^{K_1} f(1, u) du ds \\ + \frac{K_0 b_1^2 (b_2/\lambda) e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{(1 + \frac{b_2}{\lambda}(K_1 + t))(1 + \frac{b_2}{\lambda}t)} \lambda \int_0^t \int_0^{K_0} u f(0, u) du ds. \quad (4.33)$$

Hence, we can write

$$P_t f(0, x) \geq C_3(t) \int_0^{K_1} f(1, u) du + C_4(t) \int_{k_0}^{K_0} f(0, u) du, \quad (4.34)$$

where

$$C_3(t) := \lambda t \frac{b_1 e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{1 + \frac{b_2}{\lambda}(K_1 + t)}, \quad C_4(t) := \frac{K_0 b_1^2 (b_2/\lambda) e^{-t(\varrho_1 + \varrho_2(K_0 \vee K_1 + t))}}{(1 + \frac{b_2}{\lambda}(K_1 + t))(1 + \frac{b_2}{\lambda}t)} \lambda k_0 t.$$

Now, we need to normalise the measures appearing in (4.32) and (4.34) and choose $\tau$ large enough so that the resulting factors are smaller than 1 (see (4.24)). After the renormalisation, we should choose $\tau$ such that $\tau \ge 2(K_0 \vee K_1)$ and

$$C_4(\tau)(K_0 - k_0) \le 1, \quad C_3(\tau)K_1 \le 1, \quad C_2(\tau)K_1 \le 1, \quad C_1(\tau)k_0K_0(K_0 - k_0) \le 1. \quad (4.35)$$

Since the values of $K_0, K_1$ depend on $R_{\gamma}$ which is itself independent of $\tau$ (see the remark just below (4.21)), the exponential decay dominates in all the expressions. Hence, for $\tau$ large enough, there exists $c \in (0, 1]$ such that for all $(\mathfrak{e}, x) \in \mathbb{K}_{\gamma}$,

$$P_{\tau} f(\mathfrak{e}, x) \geq c \left( \int_{k_0}^{K_0} f(0, u) \frac{du}{K_0 - k_0} + \int_{0}^{K_1} f(0, u) \frac{du}{K_1} \right).$$

Hence, Assumption (A3) is satisfied for the probability measure $\nu$ defined by

$$\nu(d\mathfrak{e}, dx) := \frac{\mathbb{1}_{\{k_0 \le x \le K_0\}}}{K_0 - k_0} \delta_0(d\mathfrak{e}) dx + \frac{\mathbb{1}_{\{x \le K_1\}}}{K_1} \delta_1(d\mathfrak{e}) dx. \quad (4.36)$$
**Conclusion** Now that we have proved that Assumption **A** is satisfied, we are in the position to apply Theorem 1.1 in [7] to conclude the proof. □

**Remark 4.4.** The above computations can be used to find an explicit value for the speed of convergence $\omega$ appearing in Theorem 4.3. Indeed, [7, Prop. 3.6] gives quantitative estimates depending on the parameters $\alpha, \beta, \theta, c, d$ and $\tau$ involved in Assumptions (A1)-(A4). The resulting expression is technical, but can be derived once we have identified all the parameters involved. Indeed, for us $\beta = e^{\lambda\tau}$, $d = 1$, and $\alpha$ and $\theta$ are respectively defined in (4.21) and (4.22). Obtaining explicit values for $\tau$ and $c$ is more challenging. In fact, $\tau$ depends on the constants $k_0, K_0$ and $K_1$ (which we would also need to make explicit, see the definition of the compact set $\mathbb{K}$) and is such that (4.35) holds true. The latter involves the complicated expressions $C_1(\tau)$, $C_2(\tau)$, $C_3(\tau)$ and $C_4(\tau)$. Then, $c$ would be equal to the minimum between the four values appearing in Condition (4.35) once $\tau$ has been fixed. We chose not to develop this point further as it is lengthy and will not be needed later.

# 5 Law of large numbers

Now that the asymptotic behaviour of the mean measure has been determined, we are interested in relating this behaviour to the realisations of the process $\mathcal{Z}$. Namely, in Theorem 4.3, we obtained that $e^{-\lambda t}\mathbb{E}(\langle \mathcal{Z}_t, f \rangle)$ behaves like $\langle N, f \rangle$ as $t \to \infty$ for appropriate test functions $f$. Ideally, for a realisation of the process $\mathcal{Z}$, we would also like to obtain that $\langle \mathcal{Z}_t(\omega), f \rangle$ (once correctly renormalised) behaves like $\langle N, f \rangle$ as $t \to \infty$. Unfortunately, we are only able to show the convergence in probability stated in Theorem 1.4.

To prove Theorem 1.4, we exploit the classical martingale associated to the eigenelements $(\psi, \lambda)$ (harmonic function).

We start with a preliminary lemma divided in two parts. The first one is about deriving a bound on the first moment semigroup, which is a direct consequence of Theorem 4.3. This will be used several times in the forthcoming proofs. The second part is the $L^2$ estimate we shall use for the martingale and law of large numbers below, in the vein of many-to-two formula (or formula for forks [6]).

**Lemma 5.1.** Let $\lambda$ be the eigenvalue defined in (4.2).

i) Let $z_0 = (\mathfrak{e}_0, x_0) \in S$, $\gamma \ge 2$ and $f \in \mathcal{B}(V_\gamma)$. Then, there exists a constant $C > 0$ such that

$$M_t f(z_0) \le C e^{\lambda t} (1 + x_0^\gamma).$$

ii) There exists $C > 0$ such that for any $z_0 = (\mathfrak{e}_0, x_0) \in S$ and $g: S \to \mathbb{R}$ measurable function such that which satisfy $|g(\mathfrak{e}, x)| \le 1 + x$ for any $(\mathfrak{e}, x) \in S$, we have

$$\mathbb{E}_{\delta_{z_0}} (\langle \mathcal{Z}_t, g \rangle^2) \le C e^{\lambda t} (1 + x_0^2) + \int_0^t \int_{S^4} M_{t-s} g(z_1) M_{t-s} g(z_2) K(z, dz_1 dz_2) M_s(z_0, dz) \, ds.$$

Here $M(z_0, \cdot)$ is the measure associated to the positive semigroup $M$ and defined for any Borel measurable set $A \subset S$ by $M_s(z_0, A) := M_s 1_A(z_0)$ and

$$K(z, dz_1 dz_2) := b_1 \mathfrak{e} \left\{ 2\delta_{(0,x),(1,0)}(dz_1, dz_2) + \delta_{(1,0),(1,0)}(dz_1, dz_2) \right\} \\ + \int_0^1 b_2 x \left\{ \delta_{(0,(1-\alpha)x),(\mathfrak{e},\alpha x)}(dz_1, dz_2) + \delta_{(0,(1-\alpha)x),(1,0)}(dz_1, dz_2) + \delta_{(\mathfrak{e},\alpha x),(1,0)}(dz_1, dz_2) \right\} d\alpha.$$
Proof. (Proof of Lemma 5.1)

i) Fixing $\gamma \ge 2$, Theorem 4.3 yields for any $t \ge 0$,

$$|e^{-\lambda t} M_t f(z_0) - \psi(z_0) \langle N, f \rangle| \le C e^{-wt} V_\gamma(z_0),$$

where $w > 0$ and $C > 0$. Now, one can check that $\langle N, f \rangle < \infty$ by a direct computation since $N$ decreases exponentially. Recalling that $\psi \le V_\gamma$ and $V_\gamma(z_0)$ is dominated by $1 + x_0^\gamma$ ends the proof of i).

ii) For the $L^2$ computation, we follow [6] and use the underlying genealogy and the Ulam-Harris-Neveu notation evoked in Section 1. For $t > 0$, we have

$$\langle \mathcal{Z}_t, g \rangle = \sum_{u \in V_t} g(z_t^u),$$

where we remind that $V_t$ denotes the indexing set of the individuals alive at time $t$ and $z_t^u = (e^u, x_t^u)$ denotes the characteristics at time $t$ of individual $u \in V_t$. In addition, for $u, v \in \mathcal{U}$, $u \wedge v$ denotes the label of the most recent ancestor of $u$ and $v$, and $(vi)_{i \in \{1,2,3\}}$ are the descendants of $v$, and $u \succ v$ encodes the order in the tree (i.e. $u$ is a descendant of $v$). Having all this in mind, we first notice that

$$\langle \mathcal{Z}_t, g \rangle^2 = \sum_{u,v \in V_t} g(z_t^u)g(z_t^v) = \sum_{u \in V_t} g(z_t^u)^2 + \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} I_t(w) \quad (5.1)$$

where for any $w \in \mathcal{U}$, $b(w)$ is the time at which the individual labelled by $w$ branches (by convention it is infinite if this event does not happen) and

$$I_t(w) = \sum_{\substack{u,v \in V_t \\ i,j \in \{1,2,3\}, i \neq j \\ u \succ w_i, v \succ w_j}} g(z_t^u)g(z_t^v) = \sum_{i,j \in \{1,2,3\}, i \neq j} \left( \sum_{u \in V_t, u \succ w_i} g(z_t^u) \times \sum_{v \in V_t, v \succ w_j} g(z_t^v) \right)$$

We evaluate the expectation of each term involved in $\langle \mathcal{Z}_t, g \rangle^2$. First, in view of Lemma 5.1 with $\gamma = 2$, for any $z_0 \in S$, we have

$$\mathbb{E}_{\delta_{z_0}} \left( \sum_{u \in V_t} g(z_t^u)^2 \right) = M_t(g^2)(z_0) \le C e^{\lambda t} (1 + x_0^2). \quad (5.2)$$

Second, we deal with $\mathbb{E}_{\delta_{z_0}} (\sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} I_t(w))$. For any $w \in \mathcal{U}$ and for $i \in \{1,2,3\}$, we have

$$\mathbb{1}_{b(w)<t} \mathbb{E} \left( \sum_{u \in V_t, u \succ w_i} g(z_t^u) \mid b(w), z_{b(w)}^{wi} \right) = \mathbb{1}_{b(w)<t} M_{t-b(w)} g(z_{b(w)}^{wi}).$$

For any $w \in \mathcal{U}$, the branching property then yields

$$\mathbb{1}_{b(w)<t} \mathbb{E}_{\delta_{z_0}} (I_t(w) | \mathcal{F}_{b(w)}) = \mathbb{1}_{b(w)<t} \sum_{i \neq j \in \{1,2,3\}} M_{t-b(w)} g(z_{b(w)}^{wi}) M_{t-b(w)} g(z_{b(w)}^{wj}).$$

Combining these identities, we obtain

$$\mathbb{E}_{\delta_{z_0}} \left( \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} I_t(w) \right) = \mathbb{E}_{\delta_{z_0}} \left( \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} \sum_{i \neq j \in \{1,2,3\}} M_{t-b(w)} g(z_{b(w)}^{wi}) M_{t-b(w)} g(z_{b(w)}^{wj}) \right).$$
Recalling that $B(z)$ is the branching rate of individual of type $z$ defined in (2.3), we observe that

$$ \mathbb{E}_{\delta_{z_0}} \left( \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} f(z_{b(w-)}^w, b(w)) \right) = \int_{S \times [0,t]} B(z) f(z,s) M_s(z_0, dz) ds $$

for any $f$ non-negative measurable function defined on $S \times [0, t]$. Denoting by $p(z, dz_1 dz_2 dz_3)$ the distribution of types of the three off-springs of an individual with type $z$, we obtain

$$ \begin{aligned} & \mathbb{E}_{\delta_{z_0}} \left( \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} I_t(w) \right) \\ &= \int_0^t \int_{S^4} \sum_{i \neq j \in \{1,2,3\}} M_{t-s} g(z_i) M_{t-s} g(z_j) p(z, dz_1 dz_2 dz_3) B(z) M_s(z_0, dz) ds \end{aligned} $$

Introducing the measure $K$ defined by

$$ B(z) \int_{S^3} p(z, dz_1 dz_2 dz_3) \sum_{i \neq j \in \{1,2,3\}} f_1(z_i) f_2(z_j) = \int_{S^2} K(z, dz_1 dz_2) f_1(z_1) f_2(z_2) $$

where $f_1, f_2$ are non-negative measurable functions, we obtain

$$ \begin{aligned} & \mathbb{E}_{\delta_{z_0}} \left( \sum_{w \in \mathcal{U}} \mathbb{1}_{b(w)<t} I_t(w) \right) \\ &= \int_0^t \int_{S^3} M_{t-s} g(z_1) M_{t-s} g(z_2) K(z, dz_1 dz_2) M_s(z_0, dz) ds. \end{aligned} \quad (5.3) $$

Finally, we observe that $K$ is the kernel of binary splitting of ancestral lineages at division of the individual $z$, i.e. we choose two of the three offsprings. Its form can be explicitly derived, which ends the proof of ii) by combining (5.2) and (5.3). $\square$

We consider now the classical local martingale associated to the eigenfunction $\psi$:

$$ Y_t := \exp(-\lambda t) \langle \mathcal{Z}_t, \psi \rangle, \quad t \ge 0, $$

recalling that $\lambda$ is the eigenvalue defined in (4.2) and $\psi$ is the eigenfunction defined in (1.16). We show that it converges a.s. to a non-degenerate random variable.

**Proposition 5.2.** Under Assumption (1.6) and assuming that $\mathcal{Z}_0$ has at least one atom with probability 1, $(Y_t)_{t \ge 0}$ converges a.s, as $t \to \infty$, to a positive finite random variable $W$.

*Proof. (Proof of Proposition 5.2.)* The fact that $(Y_t)_{t \ge 0}$ is a local martingale is classical since we use here the harmonic function $\psi$. It can be directly derived from (3.2).

Now, we first prove that $(Y_t)_{t \ge 0}$ is bounded in $L^2(\Omega)$ when we start from one single individual with random and bounded initial condition $Z_0 \in S$ (Step 1). This ensures that the limit is positive with positive probability. In Step 2, we obtain that the limit is a.s. positive by standard arguments using the branching property. In Step 3, we end the proof by extending the a.s. convergence to the case of initial conditions involving several individuals and unbounded types, under Assumption 1.6.
**Step 1.** Since $\psi$ is dominated by $1+x$, we may apply Lemma 5.1-ii) to $g = \psi$. Writing $\mu_0$ the law of $Z_0$ and $\mu_0 M_s(dz) = \mathbb{E}(M_s(Z_0, dz)) = \mathbb{E}_{\delta_{Z_0}}(Z_s(dz))$ and using that $M_t\psi(z) = e^{\lambda t}\psi(z)$, it leads to

$$ \mathbb{E}_{\delta_{Z_0}} (\langle \mathcal{Z}_t, \psi \rangle^2) \le C' e^{\lambda t} (1 + x_0^2) + \int_0^t e^{2\lambda(t-s)} \int_{S^3} \psi(z_1)\psi(z_2) K(z, dz_1dz_2) \mu_0 M_s(dz) ds, $$

for some $C' \ge 0$. Adding that $K((\mathfrak{e}, x), S^2) \le cx$ and defining $f(\mathfrak{e}, x) = 1 + x^3$, we have

$$ \int_{S^3} \psi(z_1)\psi(z_2) K(z, dz_1dz_2) \mu_0 M_s(dz) \le \mu_0 M_s f. $$

Use now that $f \in \mathcal{B}(V_3)$, it comes from Lemma 5.1-i)

$$ \mathbb{E}_{\delta_{Z_0}} \langle \mathcal{Z}_t, \psi \rangle^2 \le C' e^{\lambda t} (1 + x_0^2) + C'' e^{2\lambda t} (1 + x_0^3), $$

where $C'' > 0$ and $x_0$ is the a.s. bound on the second component of $Z_0$. Starting from $\mathcal{Z}_0 = \delta_{Z_0}$, it ensures that the local martingale $Y_t$ is bounded in $L^2$ and thus converges a.s. and in $L^2$ to a finite random variable $W$. This guarantees that $\mathbb{E}(W) > 0$ and $W$ is positive with positive probability.

**Step 2.** Let us now prove that $W$ is a.s. positive using the regeneration property due to open segments. The argument given here exploits the branching property along a stopping line, in the same vein as [17]. Let us write $W^{(1,0)}$ for the limiting martingale when the initial condition is one single open segment of length $x = 0$. The previous step ensures that $\mathbb{P}(W^{(1,0)} > 0) > 0$. Besides, each open segment branches at constant rate $b_1 > 0$ through the mechanism we described as “apical branching”, giving birth to two open segments of length 0. As a consequence, restricting our attention to this subpopulation of open segments (disregarding the open segments of length zero created by “lateral branching”) we obtain a binary Yule process embedded in the original process, for which the first open segment of size zero is the root. Since the Yule process a.s. tends to infinity, for any $N \ge 1$, we can consider a finite stopping line in the original process $Z$ where we have $N$ segments of initial size 0. By the branching property, each one independently gives rise to a growth-fragmentation process, with the same law as $Z$ starting from $\delta_{(1,0)}$ and corresponding martingale limits $(W_i)_{i=1,...,N}$, independent and distributed like $W^{(1,0)}$. On the event $\{W = 0\}$, each $W_i$ has to be zero, which happens with probability $\mathbb{P}(W^{(1,0)} = 0)^N$. But the latter quantity becomes arbitrarily small as $N$ becomes large. Hence, we have $\mathbb{P}(W = 0) = 0$.

**Step 3.** Let us now consider an initial condition $\mathcal{Z}_0$ satisfying (1.6) and extend the previous result by a truncation argument. The index set of the initial individuals is given by $V_0 = \{1, ..., \langle \mathcal{Z}_0, 1 \rangle\} \subset \mathbb{N}$ and $(Z_i(0))_{i \in V_0}$ are their initial type. We introduce the branching process $\mathcal{Z}^{(i)}$ issued from the single individual $i$. Note that if this individual does not exist in the process at time 0, one can artificially choose a type for it, say (1, 1) otherwise. We also set $Y_t^{(i)} = \exp(-\lambda t) \langle \mathcal{Z}_t^i, \psi \rangle$. For any integer $k$, we introduce the event

$$ A_k = \{ \langle \mathcal{Z}_0, 1 \rangle \le k, \langle \mathcal{Z}_0, 1_{\{0,1\} \times [k, \infty)} \rangle = 0 \} $$
which allows to bound both the number of individuals and their maximal length at initial time. We consider

$$M_t^k = 1_{A_k} Y_t = 1_{A_k} \sum_{i \in V_0} Y_t^{(Z_i^k)},$$

where $Z_i^k$ is the initial condition $Z_i$ whose second component has been truncated at $k$. **Step 1** ensures that for each $i$, $Y_t^{(Z_i^k)}$ converges a.s. to a positive finite random variable. Hence, on the event $A_k$, $M_t^k$ is the sum involving a bounded and non zero number of terms $(Y_t^{(Z_i^k)})$'s and each one of them has a finite positive limit as $t \to \infty$. This ensures that $M_t^k$ converges a.s. to a finite positive limit. Adding that the sequence $1_{A_k}$ increases a.s. to 1 as $k \to \infty$ thanks to (1.6), we obtain that $Y_t$ converges a.s., as $t \to \infty$, to a positive finite random variable. This ends the proof. □

Now we are ready to prove Theorem 1.4 that describes the asymptotic empirical distribution in type and lengths in our exponential growing population.

*Proof. (Proof of Theorem 1.4.)* For convenience, we introduce the probability measure $\pi$ on $S$ defined by $\langle \pi, f \rangle := \frac{1}{2}\langle N, f \rangle$. We achieve the proof in the case of a single initial individual with bounded initial condition $Z_0 \in S$. Then, the convergence can be extended as in the previous proof.

Assume we showed that for any $g: S \to \mathbb{R}$ such that

$$\sup_{(\mathfrak{c}, x) \in \mathcal{X}} |f(\mathfrak{c}, x)|/(1+x) < \infty, \quad \langle \pi, g \rangle = 0,$$

it holds

$$\mathbb{E}_{\delta_{Z_0}} \left( e^{-2\lambda t} \langle \mathcal{Z}_t, g \rangle^2 \right) \to 0, \quad t \to \infty. \quad (5.4)$$

Then, the conclusion for any $f: S \to \mathbb{R}$ dominated by $\psi$ follows by taking $g = f - \frac{\langle \pi, f \rangle}{\langle \pi, \psi \rangle} \psi$ and by applying Proposition 5.2 to obtain the following convergence in probability:

$$\frac{\langle \mathcal{Z}_t, f \rangle}{\langle \mathcal{Z}_t, \psi \rangle} = \frac{e^{-\lambda t} \langle \mathcal{Z}_t, g \rangle}{e^{-\lambda t} \langle \mathcal{Z}_t, \psi \rangle} + \frac{\langle \pi, f \rangle}{\langle \pi, \psi \rangle} \to \frac{0}{W} + \frac{\langle \pi, f \rangle}{\langle \pi, \psi \rangle}, \quad t \to \infty.$$

It remains to check that (5.4) holds true. In view of Lemma 5.1-ii), we have

$$\mathbb{E}_{\delta_{Z_0}} \left( e^{-2\lambda t} \langle \mathcal{Z}_t, g \rangle^2 \right) \le e^{-2\lambda t} \int_0^t \int_{S^3} M_{t-s}g(z_1) M_{t-s}g(z_2) K(z, dz_1dz_2) \mu_0 M_s(dz) ds + C e^{-\lambda t},$$

where $C$ is a constant depending only on the bound of the initial condition $Z_0$. Now, for any $t > 0$, apply Theorem 4.3 for $\gamma = 2$ having in mind that $\langle \pi, g \rangle = 0$. It comes

$$|e^{-\lambda t} M_t g(z)| \le C V(z) e^{-\omega t}.$$

The above estimate together with the fact that $K(z, S^2) \le Cx$ and $V(z) \le C(1+x^2)$ yields

$$\mathbb{E} \left( e^{-2\lambda t} \langle \mathcal{Z}_t, g \rangle^2 \right) \le C \int_0^t \int_S \bar{f}(z) e^{-2\lambda s - 2\omega(t-s)} M_s(z_0, dz) ds = \int_0^t e^{-2\lambda s - 2\omega(t-s)} M_s \bar{f} ds$$

with $\bar{f}(e, x) = 1 + x^5$. Applying Lemma 5.1 with $\gamma = 5$ on $M_s \bar{f}$, we obtain that

$$\mathbb{E} \left( e^{-2\lambda t} \langle \mathcal{Z}_t, g \rangle^2 \right) \le C(1+x_0^5) \int_0^t e^{s\lambda} e^{-2\lambda s - 2\omega(t-s)} ds.$$

and (5.4) is proved. □
# 6 Parameter estimation

Our first motivation for this work was to answer the following question: Given a panorama of the state of the fungal mycelium at the end time of some growth experiment such as those carried out in [20], can we reconstruct the elongation rate $v$, and the rates of apical branching ($b_1$) and lateral branching ($b_2$) that gave rise to the observed numbers and length distributions of internal and external (or terminal) segments of filaments?

Assuming that the fusion of filaments may be neglected, we may use the model developed in this work in a direct way and calibrate it with the count and length data. Indeed, recall that to obtain Theorem 1.3, we scaled time in such a way that elongation happens at speed 1. For a general $v > 0$, the quantities $b_1$ and $b_2$ appearing in the definitions of $\lambda$ (see (1.13)), $N_1$ and $N_0$ (see (1.14-1.15)) should thus be replaced by $b_1/v$ and $b_2/v$. Combining the estimator $\hat{\Lambda}_T$ of the growth rate $\lambda$ defined in (1.18) by

$$ \hat{\Lambda}_T = \frac{\langle \mathcal{Z}_T, p_e \rangle}{\langle \mathcal{Z}_T, p_l \rangle} \approx \frac{\text{nb of filament apexes at time } T}{\text{total mycelial length at time } T}, \quad (6.1) $$

where $T > 0$ is the end time of the experiment, and the explicit formulae for the stationary length distributions of external and internal segments approximated by their empirical counterparts, a simple fitting procedure would allow us to reconstruct $v$, $b_1$ and $b_2$. Note that estimating $v$ from a temporal series of panoramas taken at sufficiently close times is rather easy and may prove to be more robust than estimating $v$ through the above procedure, since the very large number of open segments of filaments offers many realisations of the same deterministic growth process (with potential measurement errors, though).

However, it was shown in [20] that anastomosis cannot be disregarded and does have an impact on the growth properties of the network, and consequently $\lambda$ cannot be simply approximated by the ratio of the number of apexes to the total length of the mycelium appearing on the right-hand-side of (6.1). A rule of thumb led the authors of [20] to conclude that if we write $\alpha_A$ for the exponential growth rate of the observed number of apexes (i.e., open ends of filaments) and $\alpha_N$ for the exponential growth rate of the observed number of internal branch points (i.e., ends of internal segments of filaments), then these quantities can be related to the theoretical growth rate $\alpha_b = \lambda$ at which branching globally increases the total number of nodes and to the rate $\alpha_a$ at which anastomosis turns open ends of filaments into internal branch points as follows (see Equations (1) and (2) in [20]):

$$ \begin{aligned} \alpha_N &= \lambda + \alpha_a, \\ \alpha_A &= \lambda - \alpha_a. \end{aligned} \quad (6.2) $$

Therefore, using a temporal series of panoramas instead of a single “final” panorama, we may obtain a measure of $v$ as discussed above, and of $\lambda$ by writing

$$ \hat{\lambda}_T \approx \frac{\alpha_N(T) + \alpha_A(T)}{2}, \quad (6.3) $$

where $\alpha_N(T)$ (resp., $\alpha_A(T)$) are the empirical slopes of the logarithm of the number of internal branch points (resp., of apexes) through time over the time interval $[0, T]$. Once we have an empirical value for $v$ and $\lambda$, we are back in line with our simple branching model (since $\lambda$ corresponds to the growth rate of the number of nodes due to branching) and Equation (1.13) applies. But a last hurdle remains: anastomosis may have stopped the growth of some of the observed segments (and turned open segments into closed ones), distorting the stationary length profiles $N_0$ and $N_1$.
To circumvent this problem, note that in practice the centre of the mycelial network rapidly becomes very dense in such experiments, rendering the identification and measurement of segments in this region very difficult. In [20], the authors resorted to considering only the information contained in a large ring centred at the location of the initial spore but excluding the denser central part of the network (see Figure 3 in [20]). The second advantage of this approach is that, because the spatial spread of the mycelium happens radially, the directions taken by the “primary” filaments serving as a backbone for the network diverge and the network becomes more and more well-spread as we go from the centre towards the outside. This property renders the measurements easier in the ring, and in this region anastomosis mainly happens to relatively long filaments. This suggests that matching the data on types and lengths to the stationary distributions $N_1, N_0$ over an interval of lengths of the form (0, L) only, for some small L that will have to be determined empirically, may allow us to get around the fact that anastomosis prevents some of the long segments to occur and thereby distorts the tail of the distribution in segment lengths. We shall pursue this direction in future work.

**Acknowledgement.** This work is part of the NEMATIC project on "Growing and branching networks: Analysis, modelling and simulation of multiscale spatial exploration, spreading and morphogenesis under constraints". The authors thank the whole group for inspiring discussions and, in particular, Florence Chapeland-Leclerc, Eric Herbert and GWenaël Ruprich-Robert for many explanations on the physiology of filamentous fungi and very useful comments on modelling choices, Jules Depersin for some preliminary work that helped to shape the approach developed here, and Yves D'Angelo for constant efforts to keep the NEMATIC dynamics going despite the current research conditions. They warmly thank Pierre Gabriel for several useful discussions on the asymptotic behaviour of the semigroup considered in Section 4.2. Finally, they are grateful to the two anonymous reviewers for their comments which helped them to improve the presentation of the most technical parts of the paper. M.T. acknowledges financial support by the DIM Math Innov funded by the Paris region, and by the Fondation Mathématique Jacques Hadamard. V.B. was supported by the ANR project ABIM (ANR-16-CE40-0001). All authors also acknowledge the support of the chaire program Mathematical Modelling and Biodiversity (Ecole polytechnique, Muséum National d'Histoire Naturelle, Veolia Environment, Fondation X), of the CNRS "PEPS - Engineering Inspired by Nature" project BioRés and of the CNRS MITI project Dynamics of Random Expanding networks: Analysis, modelling and simulation of Multi-Scale spatial exploration, spreading and morphogenesis under constraints.

## References

[1] O. Arino. A survey of structured cell population dynamics. *Acta biotheoretica*, 43 (1-2):3-25, 1995.

[2] S. Asmussen and H. Hering. Strong limit theorems for general supercritical branching processes with applications to branching diffusions. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete*, 36(3):195-212, 1976.

[3] K. B. Athreya. Limit theorems for multitype continuous time Markov branching processes. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete*, 12(4): 320-332, 1969.
[4] W. Balmant, M. H. Sugai-Guérios, J. H. Coradin, N. Krieger, A. F. Junior, and D. A. Mitchell. A model for growth of a single fungal hypha based on well-mixed tanks in series: simulation of nutrient and vesicle transport in aerial reproductive hyphae. *PLoS One*, 10(3):e0120307, 2015.

[5] J. Banasiak, K. Pichór, and R. Rudnicki. Asynchronous exponential growth of a general structured population model. *Acta applicandae mathematicae*, 119(1):149–166, 2012.

[6] V. Bansaye, J.-F. Delmas, L. Marsalle, and V. C. Tran. Limit theorems for Markov processes indexed by continuous time Galton-Watson trees. *Ann. Appl. Probab.*, 21 (6):2263–2314, 2011.

[7] V. Bansaye, B. Cloez, P. Gabriel, and A. Marguet. A non-conservative Harris’ ergodic theorem. *To appear in J. London Mathematical Society, arXiv 1903.03946*, 2019.

[8] D. J. Barry and G. A. Williams. Microscopic characterisation of filamentous microbes: towards fully automated morphological quantification through image analysis. *Journal of microscopy*, 244(1):1–20, 2011.

[9] J. Bertoin and A. R. Watson. A probabilistic approach to spectral analysis of growth-fragmentation equations. *Journal of Functional Analysis*, 274(8):2163–2204, 2018.

[10] L. Boddy, J. Wood, E. Redman, J. Hynes, and M. D. Fricker. Fungal network responses to grazing. *Fungal Genetics and Biology*, 47(6):522–530, 2010.

[11] G. P. Boswell and F. A. Davidson. Modelling hyphal networks. *Fungal Biology Reviews*, 26(1):30–38, 2012.

[12] F. Brikci, J. Clairambault, and B. Perthame. Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle. *Mathematical and Computer Modelling*, 47(7-8):699–713, 2008.

[13] J. Calvo, M. Doumic, and B. Perthame. Long-time asymptotics for polymerization models. *Communications in Mathematical Physics*, 363(1):111–137, 2018.

[14] J. A. Cañizo, P. Gabriel, and H. Yoldas. Spectral gap for the growth-fragmentation equation via Harris’s Theorem. *SIAM Journal on Mathematical Analysis*, 53(5):5185–5214, 2021.

[15] V. Capasso and F. Flandoli. On the mean field approximation of a stochastic model of tumour-induced angiogenesis. *European Journal of Applied Mathematics*, 30(4): 619–658, 2019.

[16] R. Catellier, Y. D’Angelo, and C. Ricci. A mean-field approach to self-interacting networks, convergence and regularity. *Mathematical Models and Methods in Applied Sciences*, 31(13):2597–2641, 2021.

[17] B. Chauvin. Sur la propriété de branchement. *Ann. IHP Probab. Statist.*, 22(2): 233–236, 1986.

[18] B. Cloez. Limit theorems for some branching measure-valued processes. *Advances in Applied Probability*, 49(2):549–580, 2017.
[19] B. Cloez, B. de Saporta, and T. Roget. Long-time behavior and Darwinian optimality for an asymmetric size-structured branching process. *Journal of Mathematical Biology*, 83(6):1–30, 2021.

[20] J. Dikec, A. Olivier, C. Bobée, Y. D'Angelo, R. Catellier, P. David, F. Filaine, S. Herbert, C. Lalanne, H. Lalucque, L. Monasse, M. Rieu, G. Ruprich-Robert, A. Véber, F. Chapeland-Leclerc, and E. Herbert. Hyphal network whole field imaging allows for accurate estimation of anastomosis rates and branching dynamics of the filamentous fungus Podospora anserina. *Scientific Reports*, 10(1):3131, 2020.

[21] M. Doumic, A. Olivier, and L. Robert. Estimating the division rate from indirect measurements of single cells. *arXiv preprint*, arXiv:1907.05108, 2019.

[22] H. Du, M. Ayouz, P. Lv, and P. Perré. A lattice-based system for modeling fungal mycelial growth in complex environments. *Physica A: Statistical Mechanics and its Applications*, 511:191–206, 2018.

[23] H. Du, T.-B.-T. Tran, and P. Perré. A 3-variable PDE model for predicting fungal growth derived from microscopic mechanisms. *Journal of Theoretical Biology*, 470: 90–100, 2019.

[24] J. Engländer, S. C. Harris, and A. E. Kyprianou. Strong law of large numbers for branching diffusions. *Annales de l'I.H.P. Probabilités et statistiques*, 46(1):279–298, 2010.

[25] M. Fricker, L. Heaton, N. Jones, and L. Boddy. The mycelium as a network. *Microbiol. Spectr.*, 5(3):3131, 2017.

[26] L. Heaton, B. Obara, V. Grau, N. Jones, T. Nakagaki, L. Boddy, and M. D. Fricker. Analysis of fungal networks. *Fungal Biology Reviews*, 26(1):12–29, 2012.

[27] H. Kesten and B. P. Stigum. A limit theorem for multidimensional Galton-Watson processes. *Ann. Math. Statist.*, 37(5):1211–1223, 1966.

[28] A. Lamour, F. Van Den Bosch, A. J. Termorshuizen, and M. J. Jeger. Modelling the growth of soil-borne fungi in response to carbon and nitrogen. *IMA J. Math. Appl. Med. Biol.*, 17:329–346, 2000.

[29] A. Marguet. Uniform sampling in a structured branching population. *Bernoulli*, 25 (4A):2649–2695, 2019.

[30] S. Mischler and J. Scher. Spectral analysis of semigroups and growth-fragmentation equations. *Annales de l'Institut Henri Poincare (C) Non Linear Analysis*, 33(3):849–898, 2016.

[31] E. Nummelin. *General Irreducible Markov Chains and Non-Negative Operators*. Cambridge Tracts in Mathematics. Cambridge University Press, 1984.

[32] B. Perthame. *Transport equations in biology*. Springer Science and Business Media, 2006.

[33] S. H. Tindemans, N. Kern, and B. M. Mulder. The diffusive vesicle supply center model for tip growth in fungal hyphae. *Journal of Theoretical Biology*, 238(4):937–948, 2006.
[34] V. C. Tran. Large population limit and time behaviour of a stochastic particle model describing an age-structured population. *ESAIM Probab. Stat.*, 12:345–386, 2008.

# Appendix

## A.1 Characterisation of eigenvalues

In this section, we consider the spectral problem (4.1) with eigenfunctions $N_1$ and $N_0$ satisfying the following conditions:

$$ \int_{\mathbb{R}_+^*} (1+x)|N_0(x)|dx < \infty, \quad \int_{\mathbb{R}_+} (1+x)|N_1(x)|dx < \infty, \\ \lim_{x \to \infty} N_1(x) = 0 \text{ and } \lim_{x \to \infty} xN_1(x) = 0. $$

We prove that the associated eigenvalues satisfy (1.13).

First, observe that

$$ \int_0^\infty N_1(x)dx = \int_0^\infty N_0(x)dx. \quad (A.1) $$

To see this, it suffices to integrate (4.1a) and (4.1c) over $\mathbb{R}_+$ and then to subtract one expression from the other. Since $\lim_{x \to \infty} N_1(x) = 0$, we have

$$ (b_1+\lambda) \int_0^\infty N_1(x)dx - \lambda \int_0^\infty N_0(x)dx \\ = N_1(0) - b_1 \int_0^\infty N_1(x)dx - b_2 \int_0^\infty x(N_1(x)+N_0(x))dx. $$

Plugging in the value prescribed for $N_1(0)$ in (4.1b), we can conclude that (A.1) holds true.

Let us now sum up the two integrated equations. It comes

$$ \lim_{x \to \infty} \left( N_1(x) - N_1(0) + \lambda \int_0^\infty (N_1(x) + N_0(x))dx \right) = b_2 \int_0^\infty x(N_1(x) + N_0(x))dx. $$

Using (4.1b), we can write that

$$ \lambda \int_0^\infty (N_1(x) + N_0(x))dx = 2b_1 \int_0^\infty N_1(x)dx + 2b_2 \int_0^\infty x(N_1(x) + N_0(x))dx. \quad (A.2) $$

Multiplying both (4.1a) and (4.1c) by $x$, integrating over $\mathbb{R}_+$, summing up the two expressions just obtained and then integrating by parts and using Fubini's theorem, we arrive at

$$ [xN_1(x)]_0^\infty - \int_0^\infty N_1(x)dx + \lambda \int_0^\infty x(N_1(x)+N_0(x))dx = 0. $$

Using that $\lim_{x \to \infty} xN_1(x) = 0$, we obtain

$$ \lambda \int_0^\infty x(N_1(x) + N_0(x))dx = \int_0^\infty N_1(x)dx. \quad (A.3) $$

Plugging (A.3) into (A.2) and using (A.1), we can finally conclude that $\lambda$ satisfies the desired equation

$$ \lambda = b_1 + \frac{b_2}{\lambda}. $$
## A.2 Proof of Proposition 4.1

We first deal with (4.1a). Let us define $\tilde{N}_1(x) = (b_1 + b_2x + \lambda)e^{-\int_0^x (b_1+b_2y+\lambda)dy}$. Notice that

$$ \tilde{N}_1'(x) = (b_2 - (b_1 + b_2x + \lambda)^2)e^{-\int_0^x (b_1+b_2y+\lambda)dy} $$

and

$$ \int_x^\infty \tilde{N}_1(y)dy = e^{-\int_0^x (b_1+b_2y+\lambda)dy} $$

Hence, $\tilde{N}_1$ obviously satisfies (4.1a). Moreover, notice that $\int_0^\infty \tilde{N}_1(x)dx = 1$. Now, let us define

$$ \tilde{N}_0(x) := \frac{\tilde{N}_1(x)}{(b_2x + \lambda)^2}(b_2 + b_1(b_2x + \lambda)) + \frac{b_2e^{-\int_0^x (b_1+b_2y+\lambda)dy}}{(b_2x + \lambda)^3}(2b_2 + b_1(b_2x + \lambda)). $$

We first notice that

$$ \int_0^\infty \tilde{N}_0(x)dx = \int_0^\infty \tilde{N}_1(x)dx. \qquad (A.4) $$

Indeed, from the definition of $\tilde{N}_1$ we have that

$$ \int_0^\infty \tilde{N}_0(x)dx = \int_0^\infty \frac{\tilde{N}_1(x)}{(b_2x + \lambda)^2}(b_2 + b_1(b_2x + \lambda))dx \\ + b_2 \int_0^\infty \frac{2b_2 + b_1(b_2x + \lambda)}{(b_2x + \lambda)^3} \int_x^\infty \tilde{N}_1(y)dy dx =: I + II. $$

Integrating by parts in the second term in the above expression, we obtain that

$$ II = \left(\frac{b_1}{\lambda} + \frac{b_2}{\lambda^2}\right) \int_0^\infty \tilde{N}_0(x)dx - I. $$

Combining these identities with the fact that $\lambda$ satisfies the equality $\lambda^2 = b_1\lambda + b_2$, we obtain that (A.4) holds true. This necessarily implies that the first condition in (4.5) is satisfied.

Now, we turn to (4.1c). A tenacious reader may check that the following is true: For any $x > 0$

$$ \int_x^\infty \tilde{N}_0(y)dy = \left( \frac{b_1}{\lambda + b_2x} + \frac{b_2}{(\lambda + b_2x)^2} \right) e^{-\int_0^x (b_1 + \lambda + b_2y) dy}. $$

Hence, plugging $(\tilde{N}_0, \tilde{N}_1)$ and the above value into (4.1c) and rearranging the terms, we obtain that the couple $(\tilde{N}_0, \tilde{N}_1)$ indeed satisfies (4.1c).

In order to finish with the system (4.1), it remains to verify that the value prescribed in (4.1b) matches $\tilde{N}_1(0) = b_1 + \lambda$. Using the equation satisfied by $\lambda$ and the fact that $\int_0^\infty \tilde{N}_1(x)dx = 1$, we have

$$ \tilde{N}_1(0) = b_1 + \lambda = 2b_1 \int_0^\infty \tilde{N}_1(x)dx + \frac{b_2}{\lambda} \int_0^\infty \tilde{N}_1(x)dx. $$

To obtain the condition in (4.1b), notice that $\tilde{N}_1$ and $\tilde{N}_0$ satisfy the relationship in (A.3) thanks to the fact that they solve (4.1a) and (4.1c), respectively. Plugging this relationship in the above expression gives the desired result. We have now shown that $(N_1, N_0) = (\tilde{N}_1, \tilde{N}_0)$.
We now proceed with $\psi$. For any $c \in \mathbb{R}$, plugging

$$ \psi_1(x) = c \left( 1 + \frac{b_2}{\lambda} x \right) \quad \text{and} \quad \psi_0(x) = c \frac{b_2 x}{\lambda} $$

in (4.4a) and (4.4b), we easily obtain that these functions are solutions to this system of equations.

The relation in (4.5) fixes the above constant $c$. Indeed, we must have

$$ c \left( 1 + \frac{b_2}{\lambda} \int x(N_1(x) + N_0(y)) dx \right) = 1. $$

From (A.3), we obtain that $c = \frac{\lambda^2}{\lambda^2 + b_2}$.