# A Mean-field forest-fire model

Xavier Bressaud, Nicolas Fournier

► To cite this version:

Xavier Bressaud, Nicolas Fournier. A Mean-field forest-fire model. ALEA : Latin American Journal of Probability and Mathematical Statistics, 2014, 11 (2), pp.589-614. hal-01307369

HAL Id: hal-01307369

https://hal.sorbonne-universite.fr/hal-01307369

Submitted on 26 Apr 2016

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
# A Mean-field forest-fire model

Xavier Bressaud and Nicolas Fournier

Université Paul Sabatier, Institut de Mathématiques de Toulouse, F-31062 Toulouse Cedex,
France

E-mail address: bressaud@math.univ-toulouse.fr

URL: http://www.math.univ-toulouse.fr/~bressaud/

Laboratoire de Probabilités et Modèles Aléatoires, CNRS UMR 7599, Université Paris 6,
Cas courrier 188, 4 place Jussieu, 75252 Paris Cedex 05, France.

E-mail address: nicolas.fournier@upmc.fr

URL: http://www.proba.jussieu.fr/pageperso/fournier/

**Abstract.** We consider a family of discrete coagulation-fragmentation equations closely related to the one-dimensional forest-fire model of statistical mechanics: each pair of particles with masses $i, j \in \mathbb{N}$ merge together at rate 2 to produce a single particle with mass $i+j$, and each particle with mass $i$ breaks into $i$ particles with mass 1 at rate $(i-1)/n$. The (large) parameter $n$ controls the rate of ignition and there is also an acceleration factor (depending on the total number of particles) in front of the coagulation term. We prove that for each $n \in \mathbb{N}$, such a model has a unique equilibrium state and study in details the asymptotics of this equilibrium as $n \to \infty$: (I) the distribution of the mass of a typical particle goes to the law of the number of leaves of a critical binary Galton-Watson tree, (II) the distribution of the mass of a typical size-biased particle converges, after rescaling, to a limit profile, which we write explicitly in terms of the zeroes of the Airy function and its derivative. We also indicate how to simulate perfectly a typical particle and a size-biased typical particle by pruning some random trees.

## 1. Introduction

1.1. *The forest-fire model.* The forest-fire model of statistical mechanics has been introduced by Henley (1989) and Drossel and Schwabl (1992) in the context of self-organized criticality. From the rigorous point of view, the one-dimensional forest-fire model has been studied by van den Berg and Járai (2005), Brouwer and Pennanen (2006) and by the authors Bressaud and Fournier (2010, 2013). We refer to the introduction of Bressaud and Fournier (2013) for many details.

<footnote>Received by the editors June 20, 2014; accepted October 27, 2014.</footnote>

<footnote n="2010">2010 Mathematics Subject Classification. 82B05.</footnote>

<footnote>Key words and phrases. Self-organized criticality, Smoluchowski's equation, Coagulation, Fragmentation, Equilibrium Asymptotic Behavior, Forest-fire model, Galton-Watson trees, Continuum random trees, Pruning, Scaling limit.</footnote>
Let us now describe the one-dimensional forest-fire model: on each site of $\mathbb{Z}$, seeds fall at rate 1 and matches fall at rate $1/n$, for some (large) $n \in \mathbb{N}$. Each time a seed falls on a vacant site, this site immediately becomes occupied (by a tree). Each time a match falls on an occupied site, it immediately burns the corresponding occupied connected component.

From the point of view of self-organized criticality, it is interesting to study what happens when $n$ increases to infinity. Then matches are very rare, but tree clusters are huge before they burn. In Bressaud and Fournier (2010), we have established that after normalization, the forest-fire process converges, as $n \to \infty$, to a scaling limit.

## 1.2. A related mean-field model.
We now introduce a mean-field model formally related to the forest-fire process, see Bressaud and Fournier (2009, Section 6) for a similar study when $n = 1$. Assume that each edge of $\mathbb{Z}$ has mass 1. Say that two adjacent edges $(i-1, i)$ and $(i, i+1)$ are glued if the site $i$ is occupied. Then adjacent clusters coalesce at rate 1 and each cluster with mass $k$ (containing $k$ edges and $k-1$ sites) breaks up into $k$ clusters with mass 1 at rate $(k-1)/n$.

We assign mass to edges rather than sites to preserve mass at each event. Assigning mass to sites would lead to non conservative dynamics.

Denote by $c_k^n(t)$ the concentration (number per unit of length) of clusters with mass $k \ge 1$ at time $t \ge 0$. Then the total mass should satisfy $\sum_{k \ge 1} kc_k^n(t) = 1$ for all $t \ge 0$. Neglecting correlations (which is far from being justified), the family $(c_k^n(t))_{t \ge 0, k \ge 1}$ would satisfy the following system of differential equations called (CF$_n$):

$$ \frac{d}{dt}c_1^n(t) = -2c_1^n(t) + \frac{1}{n}\sum_{k \ge 2} k(k-1)c_k^n(t), \qquad (1.1) $$

$$ \frac{d}{dt}c_k^n(t) = -(2+(k-1)/n)c_k^n(t) + \frac{1}{\sum_{l \ge 1} c_l^n(t)} \sum_{i=1}^{k-1} c_i^n(t)c_{k-i}^n(t) \quad (k \ge 2). \quad (1.2) $$

The first term on the RHS of (1.1) expresses that a cluster with mass 1 disappears at rate 2 (because it glues with each of its two neighbors at rate 1); the second term on the RHS of (1.1) says that $k$ clusters with mass 1 appear each time a cluster with mass $k$ takes fire, which occurs at rate $(k-1)/n$. The first term on the RHS of (1.2) explains that a cluster with mass $k$ disappears at rate $2+(k-1)/n$ (because it glues with each of its two neighbors at rate 1 and takes fire at rate $(k-1)/n$). Finally, the second term on the RHS of (1.2) says that when a seed falls between two clusters with masses $i$ and $k-i$, a cluster with mass $k$ appears. The number per unit of length of pairs of neighbor clusters with masses $i$ and $k-i$ is nothing but $c_i^n(t)c_{k-i}^n(t)/\sum_{l \ge 1} c_l^n(t)$. Here we implicitly use an independence argument which is not valid for the true forest-fire model.

The system (1.1)-(1.2) can almost be seen as a special coagulation-fragmentation equation, see e.g. Aizenman and Bak (1979) and Carr (1992), where particles with masses $i$ and $j$ coalesce at constant rate $K(i, j) = 2$ and where particles with mass $i \ge 2$ break up into $i$ particles with mass 1 at rate $F(i; 1, \dots, 1) = (i-1)/n$. However there is the acceleration factor $1/\sum_{l \ge 1} c_l^n(t)$ in front of the coagulation term.
1.3. On the link between the two models. The link between ($CF_n$) and the forest-fire model is only formal. Observe however that if there are only seeds **or** only matches, then the link is rigorous.

(i) Assume that all sites of $\mathbb{Z}$ are initially vacant and that seeds fall on each site of $\mathbb{Z}$ at rate 1, so that each site is occupied at time $t$ with probability $1 - e^{-t}$. Call $p_k(t)$ the probability that the edge (0, 1) belongs to a cluster with mass $k$ at time $t$. A simple computation shows that $p_k(t) = k(1 - e^{-t})^{k-1}e^{-2t}$. By space stationarity, the concentration $c_k(t)$ of particles with mass $k$ at time $t$ satisfies $c_k(t) = p_k(t)/k = (1 - e^{-t})^{k-1}e^{-2t}$. Then one easily checks that the family $(c_k(t))_{k \ge 1, t \ge 0}$ satisfies (1.1)-(1.2) with no fragmentation term (i.e. $n = \infty$).

(ii) The fragmentation term is linear and generates a priori no correlation. Assume that the successive masses of the clusters are initially distributed as a stationary renewal process with concentrations $(c_k(0))_{k \ge 1}$. In particular, the edge (0, 1) belongs to a cluster with mass $k$ with probability $kc_k(0)$. Assume also that only matches fall, at rate $1/n$ on each site. Then the probability $p_k(t)$ that the edge (0, 1) belongs to a cluster with mass $k$ at time $t$ is simply given by $p_k(t) = kc_k(0)e^{-(k-1)t/n}$ if $k \ge 2$ (here $e^{-(k-1)t/n}$ is the probability that no match has fallen on our cluster before time $t$) and $p_1(t) = c_1(0) + \sum_{k \ge 2} kc_k(0)[1 - e^{-(k-1)t/n}]$. Writing $c_k(t) = p_k(t)/k$ as previously, we see that the family $(c_k(t))_{k \ge 1, t \ge 0}$ satisfies the fragmentation equations $\frac{d}{dt}c_1(t) = n^{-1} \sum_{k \ge 2} k(k-1)c_k(t)$ and, for $k \ge 2$, $\frac{d}{dt}c_k(t) = -n^{-1}(k-1)c_k(t)$.

When one takes into account both coalescence and fragmentation, the rigorous link between the two models breaks down: in the true forest-fire model, fragmentation (fires) produces small clusters which are close to each other, so that a small cluster has more chance to have small clusters as coalescence partners. However, we have seen numerically in Bressaud and Fournier (2009, Section 6) that at equilibrium, in the special case where $n = 1$, the two models are very close to each other.

1.4. *Motivation.* Initially, the motivation of the present study was to decide if (1.1)-(1.2) is a good approximation of the true forest-fire model, at least from a qualitative point of view: do we have the same scales and same features (as $n \to \infty$)? We will see that this is not really the case. However, we believe that (1.1)-(1.2) is a very interesting model, at least theoretically, since

• many explicit computations are possible for (1.1)-(1.2), as $n \to \infty$,

• we observe self-organized criticality,

• we show that two interesting points of view (size-biased and non size-biased particles' mass distribution) lead to quite different conclusions.

1.5. *Summary of the main results of the paper.* We show in this paper the existence of a unique equilibrium state $(c_k^n)_{k \ge 1}$ with total mass $\sum_{k \ge 1} kc_k^n = 1$ for ($CF_n$), for each $n \ge 1$ fixed and we study the asymptotics of rare fires $n \to \infty$.

(I) We show that the particles' mass distribution ($p_k^n$)$_{k \ge 1}$, defined by $p_k^n = c_k^n / \sum_{l \ge 1} c_l^n$, goes weakly to the law ($p_k$)$_{k \ge 1}$ of the number of leaves of a critical binary Galton-Watson tree, which is explicit and satisfies $p_k \sim (2\sqrt{\pi}k^{3/2})^{-1}$ as $k \to \infty$;
(II) We prove that the size-biased particles’ mass distribution $(kc_k^n)_{k \ge 1}$ goes weakly to, after normalization of the masses by $n^{-2/3}$, to a continuous limit profile $(xc(x))_{x \in (0, \infty)}$ with total mass 1, for which we have two explicit expressions: the first one involves the zeroes of the Airy function and its derivative, while the Laplace transform of the Brownian excursion’s area appears in the second one. Furthermore, we check that $c(x) \sim \kappa_1 x^{-3/2}$ as $x \to 0$ and $c(x) \sim \kappa_2 e^{-\kappa_3 x}$ as $x \to \infty$, for some positive explicit constants $\kappa_1, \kappa_2, \kappa_3$.

(III) We also explain how to simulate perfectly, for $n \ge 1$ fixed, a random variable $X_n$ with law $(p_k^n)_{k \ge 1}$ and a random variable $Y_n$ with law $(kc_k^n)_{k \ge 1}$, using some pruned Galton-Watson trees.

Let us discuss briefly points (I) and (II). The particles’ mass distribution is, roughly, the law of the mass of a particle chosen uniformly at random. The size-biased particles’ mass distribution is, roughly, the law of the mass of the particle containing a given *atom*, this atom being chosen uniformly at random (think that a particle with mass *k* is composed of *k* atoms).

Point (I) says that if one picks a particle at random, then its mass $X_n$ is finite (uniformly in $n$) and goes in law, as $n \to \infty$, to a *critical* probability distribution (with infinite expectation). Point (II) says that if one picks an atom at random, then the mass $Y_n$ of the particle including this atom is of order $n^{2/3}$ and $n^{-2/3}Y_n$ goes in law to an explicit probability distribution with moments of all orders.

1.6. *Comments*. Let us now comment on these results.

(a) Self-organized criticality, see Bak-Tang-Wiesenfield Bak et al. (1987) and Henley (1989), is a popular concept in physics. The main idea is the following: in statistical mechanics, there are often some *critical* parameters, for which special features occur. Consider e.g. the case of percolation in $\mathbb{Z}^2$, see Grimmett (1989): for $p \in (0,1)$ fixed, open each edge independently with probability $p$. There is a critical parameter $p_c \in (0,1)$ such that for $p \le p_c$, there is a.s. no infinite open path, while for $p > p_c$, there is a.s. one infinite open path. If $p < p_c$, all the open paths are small (the cluster-size distribution has an exponential decay). If $p > p_c$, there is only one infinite open cluster, which is huge, and all the finite open paths are small. But if $p = p_c$, there are some large finite open paths (with a heavy tail distribution). And it seems that such phenomena, reminiscent from criticality, sometimes occur in nature, where nobody is here to finely tune the parameters. Hence one looks for models in which criticality occurs naturally. We refer to the introduction of Bressaud and Fournier (2013) for many details.

(b) Thus we observe here self-organized criticality for the particles’ mass distribution, since the limit distribution $(p_k)_{k \ge 1}$ has a heavy tail and is indeed related to a *critical* binary Galton-Watson process. Observe that point (I) is quite strange at first glance. Indeed, when $n \to \infty$, the time-dependent equations (1.1)-(1.2) tend to some coagulation equations without fragmentation, for which there is no equilibrium (because all the particles’ masses tend to infinity). However, the equilibrium state for (1.1)-(1.2) tends to some non-trivial equilibrium state as $n \to \infty$. This is quite surprising: the limit (as $n \to \infty$) of the equilibrium is not the equilibrium of the limit. Observe that this is indeed *self-organized criticality*: the *critical* Galton-Watson tree appears automatically in the limit $n \to \infty$. A possible heuristic argument is that in some sense, in the limit $n \to \infty$, only infinite clusters are
destroyed. We thus let clusters grow as much as they want, but we do not let them become infinite. Thus the system reaches *by itself* a critical state.

(c) For the size-biased particles' mass distribution, we observe no self-organized criticality, since the limit profile has an exponential decay. The two points of view (size-biased or not) seem interesting and our results show that they really enjoy different features. Let us insist on the fact that the particles' mass distribution converges without rescaling, while the size-biased particles' mass distribution converges after rescaling. This is due to the fact that there are many small particles and very few very large particles, so that when one picks a particle at random, we get a rather small particle, while when one picks an atom at random, it belongs to a rather large particle. Mathematically, since the limiting particles' mass distribution has no expectation, it is no more possible to write properly the corresponding size-biased distribution: a normalization is necessary.

(d) Let us mention the paper of Ráth and Tóth (2009), who consider a forest-fire model on the complete graph. In a suitable regime, they obtain the same critical distribution ($p_k”_{k≥1}$ as we do (see Ráth and Tóth (2009, Formula (14))). But in their case, this is the limit of the *size-biased particles' mass distribution*. Their model rather corresponds to the case of a multiplicative coagulation kernel (clusters of masses $k$ and $l$ coalesce at rate $kl$). Hence we exhibit, in some sense, a link between the Smoluchowski equation with constant kernel and the Smoluchowski equation with multiplicative kernel. See Remark 3.4 for a precise statement. In the same spirit, recall the link found by Deaconu and Tanré (2000) between multiplicative and additive coalescence.

(e) For the true forest-fire process, we proved in Bressaud and Fournier (2010) the presence of *macroscopic* clusters, with masses of order $n/\log n$ and of *microscopic* clusters, with masses of order $n^z$, for all values of $z \in [0, 1)$. Here the scales are thus very different. But there might be a (quite unclear) similarity: the asymptotic size-biased particle's mass distribution is singular at 0 (which expresses the presence of very small particles).

(f) The trend to equilibrium for coagulation-fragmentation equations has been much studied, see e.g. Aizenman and Bak (1979), Whittle (1986) and Carr (1992), under a reversibility condition, often called *detailed balance condition*. Such a reversibility assumption cannot hold here, because particles merge by pairs and break into an arbitrary large number of smaller clusters. Without reversibility, much less is known: a special case has been studied by Dubovskiǐ and Stewart (1996) and a general result has been obtained in Fournier and Mischler (2004) under a smallness condition saying that fragmentation is much stronger than coalescence. None of the above results may apply to the present model (at least for $n$ slightly large). For the specific model under study we are only able to prove the existence and uniqueness of the equilibrium. Despite much effort, we have not been able to check the convergence to equilibrium.

1.7. *Outline of the paper.* The next section is devoted to the precise statement of points (I) and (II). Section 3 contains the proofs of points (I) and (II), which are purely analytic. Using discrete pruned random trees, we indicate how to simulate perfectly a typical particle and a size-biased typical particle in Section 4.
1.8. Probabilistic interpretation of point (II). Let us finally mention that in an unpublished longer version of the present paper Bressaud and Fournier (2014), we tried to understand why the Brownian excursion arises in point (II), see Theorem 2.4 below. The main idea is that that the pruned Galton-Watson tree of which the number of leaves is $(kc_k^n)_{k \ge 1}$-distributed has a scaling limit, which is nothing but a pruned version of the famous self-similar continuum random tree of Aldous (1991a,b, 1993). The CRT is closely linked to the Brownian excursion.

Our pruning procedure is as follows: we consider a self-similar CRT, we choose leaves at random according to a Poisson measure with intensity $d(z)\nu(dz)$, where $\nu$ is the uniform measure on leaves and, for $z$ a leave, $d(z)$ is its distance from the root. These leaves send a cut-point, chosen uniformly on its branch (joining it to the root). Then we prune according to these cut-points, in a suitable order. The cut-points sent by leaves belonging to subtrees previously pruned are deactivated.

We have observed two noticeable features, although we do not really know what to do with them. (i) Conditionally on the pruned CRT, the final cut-points are Poisson distributed on the boundary. (ii) The contour of the pruned CRT is a diffusion process, which is quite noticeable and relies on our special pruning procedure. Furthermore, this diffusion process enjoys the strange property that its drift coefficient equals the Laplace exponent of its inverse local time at 0.

## 2. Precise statements of the results

For a $[0, \infty)$-valued sequence $u = (u_k)_{k \ge 1}$ and for $i \ge 0$, we put

$$m_i(u) := \sum_{k \ge 1} k^i u_k.$$

**Definition 2.1.** Let $n \in \mathbb{N}$. A sequence $c^n = (c_k^n)_{k \ge 1}$ of nonnegative real numbers is said to solve $(E_n)$ if it is an equilibrium state for $(CF_n)$ with total mass 1:

$$m_1(c^n) = 1, \qquad (2.1)$$

$$\left(2 + \frac{k-1}{n}\right) c_k^n = \frac{1}{m_0(c^n)} \sum_{i=1}^{k-1} c_i^n c_{k-i}^n \quad (k \ge 2). \qquad (2.2)$$

Observe that it then automatically holds that

$$2c_1^n = \frac{1}{n} \sum_{k \ge 2} k(k-1)c_k^n = \frac{m_2(c^n) - 1}{n}, \qquad (2.3)$$

which is the stationary version of (1.1).

To check this last claim, multiply (2.2) by $k$, sum for $k \ge 2$, use (2.1) and that $\sum_{k \ge 2} k \sum_{i=1}^{k-1} c_i^n c_{k-i}^n = \sum_{k \ge 2} \sum_{i=1}^{k-1} (i + k - i) c_i^n c_{k-i}^n = 2 \sum_{k,l \ge 1} k c_k^n c_l^n = 2m_0(c^n)$. One finds $2(1 - c_1^n) + [m_2(c^n) - 1]/n = 2$, from which (2.3) readily follows.

To state our main results, we need some background on the Airy function Ai. We refer to Janson (2007, p 94) and the references therein. Recall that for $x \in \mathbb{R}$, $\mathrm{Ai}(x) = \pi^{-1} \int_{0}^{\infty} \cos(t^3/3 + xt)dt$ is the unique solution, up to normalization, to the differential equation $\mathrm{Ai}''(x) = x\mathrm{Ai}(x)$ that is bounded for $x \ge 0$. It extends to an entire function. All the zeroes of the Airy function and its derivative lie on negative real axis. Let us denote by $a'_1 < 0$ the largest negative zero of $\mathrm{Ai}'$ and by
· · · < a$_{3}$ < a$_{2}$ < a$_{1}$ < 0 the ordered zeroes of Ai. We know that | a$_{1}$ | ≃ 2 . 338 and | a ′ $_{1}$| ≃ 1 . 019, see Finch (2004). We also know that | a$_{j}$ | ∼ (3 πj/ 2) 2 / 3 as j → ∞ , see Janson (2007, p 94).

**Theorem 2.2.** (i) For each $n \in \mathbb{N}$, $(E_n)$ has a unique solution $c^n = (c_k^n)_{k \ge 1}$.
(ii) As $n \to \infty$, there hold

$$
m_0(c^n) \sim \frac{1}{|a'_1|n^{1/3}}, \quad m_2(c^n) \sim \frac{n^{2/3}}{|a'_1|}.
$$

In some sense, $m_0(c^n)$ stands for the total concentration and $m_2(c^n)$ stands for the mean mass of (size-biased) clusters. For any fixed $l \ge 0$, we also have shown that $m_{l+1}(c^n) \sim M_l n^{2l/3}$ for some positive constant $M_l$, see Lemma 3.5 below.

Let us explain roughly these scales. Since $m_2(c^n)$ is the mean of the mass of a typical (size-biased) particle while $m_3(c^n)$ is the mean of its square, one might expect that $m_3(c^n) \simeq (m_2(c^n))^2$. By $\simeq$, we mean nothing rigorous, only that $m_3(c^n)$ and $(m_2(c^n))^2$ might have the same order of magnitude as $n \to \infty$. But some easy computations using (2.1), (2.2) and (2.3) show that $m_2(c^n) = (n-1)m_0(c^n)$ and $m_3(c^n) = 2m_2(c^n) - 1 + 2n/m_0(c^n)$. Assuming that $m_0(c^n)$ is small and that $m_2(c^n)$ is large, we thus find $m_2(c^n) \simeq nm_0(c^n)$ and $(m_2(c^n))^2 \simeq m_3(c^n) \simeq n/m_0(c^n)$, from which one easily concludes that $m_2(c^n) \simeq n^{2/3}$ and $m_0(c^n) \simeq n^{-1/3}$.

**Theorem 2.3.** For each $n \ge 1$, consider the unique solution $(c_k^n)_{k \ge 1}$ to $(E_n)$ and the corresponding particles' mass probability distribution $(p_k^n)_{k \ge 1}$ defined by $p_k^n = c_k^n/m_0(c^n)$. There holds

$$
\lim_{n \to \infty} \sum_{k \ge 1} |p_k^n - p_k| = 0, \quad \text{where} \quad p_k := \frac{2}{4^k k} \binom{2k-2}{k-1}.
$$

The sequence $(p_k)_{k \ge 1}$ is the unique nonnegative solution to

$$
\sum_{k \ge 1} p_k = 1, \quad p_k = \frac{1}{2} \sum_{i=1}^{k-1} p_i p_{k-i} \quad (k \ge 2). \tag{2.4}
$$

There holds, as $k \to \infty$,

$$
p_k \sim \frac{1}{2\sqrt{\pi k^{3/2}}}. \tag{2.5}
$$

Formally, divide (2.2) by $m_0(c^n)$ and make $n$ tend to infinity: one gets $2p_k = \sum_{i=1}^{k-1} p_i p_{k-i}$ for all $k \ge 2$. What is much more difficult (and quite surprising) is to establish that no mass is lost at the limit. Observe that (2.4) may be rewritten $\sum_{i \ge 1} p_i p_k = (1/2) \sum_{i=1}^{k-1} p_i p_{k-i}$ (for all $k \ge 2$), which corresponds to an equilibrium for a coagulation equation with constant kernel. This is quite strange, since coagulation is a monotonic process, for which no equilibrium should exist. The point is that in some sense, infinite particles are broken into particles with mass 1, in such a way that $\sum_{k \ge 1} p_k = 1$. Finally, we mention that $(p_k)_{k \ge 1}$ is the law of the number of leaves of a critical binary Galton-Watson tree, which will be interpreted in Section 4.
**Theorem 2.4.** For each $n \ge 1$, consider the unique solution $(c_k^n)_{k \ge 1}$ to $(E_n)$ and the corresponding size-biased particles' mass probability distribution $(kc_k^n)_{k \ge 1}$. For any $\phi \in C([0, \infty))$ with at most polynomial growth, there holds

$$ \lim_{n \to \infty} \sum_{k \ge 1} \phi(n^{-2/3}k)kc_k^n = \int_0^\infty \phi(x)xc(x)dx, $$

where the profile $c: (0, \infty) \mapsto (0, \infty)$ is defined, for $x > 0$, by

$$ c(x) = |a'_1|^{-1} \exp(|a'_1|x) \sum_{j=1}^{\infty} \exp(-|a_j|x). $$

The profile $c$ is of class $C^\infty$ on $(0, \infty)$, has total mass $\int_0^\infty xc(x)dx = 1$ and

$$ c(x) \stackrel{x \to 0}{\sim} \frac{1}{2\sqrt{\pi}|a'_1|x^{3/2}} \quad \text{and} \quad c(x) \stackrel{x \to \infty}{\sim} |a'_1|^{-1} \exp((|a'_1| - |a_1|)x). $$

For any $\phi \in C^1([0, \infty))$ such that $\phi$ and $\phi'$ have at most polynomial growth,

$$ 2|a'_1| \int_0^\infty \int_0^\infty x[\phi(x+y) - \phi(x)]c(x)c(y)dydx = \int_0^\infty x^2[\phi(x) - \phi(0)]c(x)dx. \quad (2.6) $$

Denote by $B_{ex}$ is the integral of the normalized Brownian excursion, see Revuz and Yor (1999, Chapter XII). For all $x > 0$,

$$ c(x) = \frac{\exp(|a'_1|x)}{2\sqrt{\pi}|a'_1|x^{3/2}} \mathbb{E} \left[ e^{-\sqrt{2}x^{3/2}B_{ex}} \right]. $$

Finally, for all $q \in (a_1 - a'_1, \infty)$ (recall that $a_1 - a'_1 < 0$),

$$ \ell(q) := \int_{0}^{\infty} (1 - e^{-qx}) c(x) dx = \frac{-\mathrm{Ai}'(q + a'_{1})}{|a'_{1}| \mathrm{Ai}(q + a'_{1})}. \quad (2.7) $$

Since the mean mass of a typical (size-biased) particle is of order $n^{2/3}$ by Theorem 2.2-(ii), it is natural to rescale the particles' masses by a factor $n^{-2/3}$. Here we state that under this scale, there is indeed a limit profile and we give some information about this profile.

### 3. Proofs

For each $n \in \mathbb{N}$, we introduce the sequence $(\alpha_k^n)_{k \ge 1}$, defined recursively by

$$ \alpha_1^n = 1, \quad (2 + (k-1)/n) \alpha_k^n = \sum_{i=1}^{k-1} \alpha_i^n \alpha_{k-i}^n \quad (k \ge 2). \quad (3.1) $$

We also introduce its generating function $f_n$, defined for $q \ge 0$ by

$$ f_n(q) = \sum_{k \ge 1} \alpha_k^n q^k. \quad (3.2) $$

Obviously, $f_n$ is increasing on $[0, \infty)$ and takes its values in $[0, \infty) \cup \{\infty\}$. The unique solution to $(E_n)$ can be expressed in terms of this sequence.
**Lemma 3.1.** Fix $n \in \mathbb{N}$. Assume that there exists a (necessarily unique) $q_n > 0$ such that $f_n(q_n) = 1$. Assume furthermore that $f'_n(q_n) < \infty$. Then there is a unique solution to $(E_n)$ and it is given by

$$c_k^n = \alpha_k^n q_n^{k-1} / f'_n(q_n) \quad (k \ge 1). \qquad (3.3)$$

Furthermore, there holds $m_0(c^n) = 1/(q_n f'_n(q_n))$.

*Proof:* We break the proof into 3 steps.

**Step 1.** A simple computation shows that for any fixed $r > 0$, $x > 0$, the sequence defined recursively by

$$y_1 = x, \quad (2 + (k - 1)/n) y_k = \frac{1}{r} \sum_{i=1}^{k-1} y_i y_{k-i} \quad (k \ge 2)$$

is given by $y_k = \alpha_k^n x(x/r)^{k-1}$.

**Step 2.** Consider a solution $(c_k^n)_{k \ge 1}$ to $(E_n)$. Then due to (2.2) and Step 1 (write $r = m_0(c^n)$ and $x = c_1^n$), $c_k^n = \alpha_k^n c_1^n (c_1^n/m_0(c^n))^{k-1}$. We deduce that

$$1 = \frac{1}{m_0(c^n)} \sum_{k \ge 1} c_k^n = \sum_{k \ge 1} \alpha_k^n (c_1^n / m_0(c^n))^k = f_n(c_1^n / m_0(c^n)).$$

Consequently, $c_1^n/m_0(c^n) = q_n$, whence $c_k^n = \alpha_k^n c_1^n q_n^{k-1}$. Next we know that $m_1(c^n) = 1$, so that $f'_n(q_n) = \sum_{k \ge 1} k \alpha_k^n q_n^{k-1} = m_1(c^n)/c_1^n = 1/c_1^n$. Consequently, $c_1^n = 1/f'_n(q_n)$ and thus $c_k^n = \alpha_k^n q_n^{k-1}/f'_n(q_n)$ as desired.

**Step 3.** Let us finally check that $c^n$ as defined by (3.3) is indeed solution to $(E_n)$ and that it satisfies $m_0(c^n) = 1/(q_n f'_n(q_n))$. First,

$$m_0(c^n) = \frac{1}{f'_n(q_n)} \sum_{k \ge 1} \alpha_k^n q_n^{k-1} = \frac{f_n(q_n)}{q_n f'_n(q_n)} = \frac{1}{q_n f'_n(q_n)}.$$

Next, (2.1) holds, since

$$m_1(c^n) = \frac{1}{f'_n(q_n)} \sum_{k \ge 1} k \alpha_k^n q_n^{k-1} = \frac{f'_n(q_n)}{f'_n(q_n)} = 1.$$

Rewriting $c_k^n = \alpha_k^n x(x/r)^{k-1}$ with $x = 1/f'_n(q_n)$ and $r = 1/(q_n f'_n(q_n))$, Step 1 implies that for $k \ge 2$,

$$(2 + (k - 1)/n) c_k^n = q_n f'_n(q_n) \sum_{i=1}^{k-1} c_i^n c_{k-i}^n = \frac{1}{m_0(c^n)} \sum_{i=1}^{k-1} c_i^n c_{k-i}^n,$$

whence (2.2). □

To go on, we need some background on Bessel functions of the first kind. Recall that for $k \ge 0$ and $z \in \mathbb{C}$,

$$J_k(z) = \frac{z^k}{2^k} \sum_{l \ge 0} \frac{(-1)^l z^{2l}}{4^l l! (l+k)!} \qquad (3.4)$$
We have the following recurrence relations, see Andrews et al. (1999, Section 4.6): for $k \ge 0$, $x \in \mathbb{R}$,

$$J'_{k}(x) = \frac{k}{x} J_{k}(x) - J_{k+1}(x) \qquad (3.5)$$

$$J_{k+2}(x) = \frac{2(k+1)}{x} J_{k+1}(x) - J_{k}(x), \qquad (3.6)$$

$$\frac{d}{dx}(x^{k+1} J_{k+1}(x)) = x^{k+1} J_{k}(x). \qquad (3.7)$$

It is known, see Andrews et al. (1999, Section 4.14), that all the zeroes of $J_k$ and $J'_k$ are real. For all $k \ge 1$, we denote by $j_k$ the first positive zero of $J_k$ and by $j'_k$ the first positive zero of $J'_k$. The sequence $(j_k)_{k \ge 0}$ is increasing, see Andrews et al. (1999, Section 4.14). Furthermore, $0 < j'_k < j_k$ for all $k \ge 1$, see Finch (2003, page 3). We will also use that there exists a constant $C > 0$ such that for all $k \ge 1$, see Finch (2003, pages 2-3),

$$k + 2^{-1/3} |a'_{1}| k^{1/3} < j'_{k} < k + 2^{-1/3} |a'_{1}| k^{1/3} + Ck^{-1/3}, \qquad (3.8)$$

where $a'_1$ is, as previously defined, the largest negative zero of $Ai'$.

**Lemma 3.2.** Let $n \in \mathbb{N}$ be fixed. The radius of convergence of the power series $f_n$ defined by (3.2) is $r_n = (j_{2n-1}/(2n))^2/2$. For all $x \in [0, r_n)$, there holds

$$f_n(x) = \sqrt{2x} \frac{J_{2n}(2n\sqrt{2x})}{J_{2n-1}(2n\sqrt{2x})}. \qquad (3.9)$$

There exists a unique $q_n \in (0, r_n)$ such that $f_n(q_n) = 1$. There holds $f'_n(q_n) = (n/q_n)(2q_n - 1 + 1/n) \in (0, \infty)$. Finally, as $n \to \infty$,

$$q_n = \frac{1}{2} [1 + |a'_1| n^{-2/3}] + O(n^{-1}).$$

*Proof:* We fix $n \in \mathbb{N}$ and divide the proof into five steps.

**Step 1.** Put $s_n = (j_{2n-1}/(2n))^2/2$. By definition of $j_{2n-1}$ and since $J_{2n-1}$ is odd on $\mathbb{R}$ and has no complex zeroes, $J_{2n-1}(2n\sqrt{2z})$ does not vanish on $\{0 < |z| < s_n\} \subset \mathbb{C}$. We thus may define $g_n(z) := \sqrt{2z}J_{2n}(2n\sqrt{2z})/J_{2n-1}(2n\sqrt{2z})$ on $\{0 < |z| < s_n\}$. Using (3.4), one easily checks that, as $z \to 0$, $g_n(z) \sim z$, so that finally, $g_n(z)$ is holomorphic on the disc $\{|z| < s_n\}$. Write $g_n(z) = \sum_{k \ge 0} \beta_k^n z^k$. Since $g_n(z) \sim z$ near 0, we deduce that $\beta_0^n = 0$ and $\beta_1^n = 1$.

**Step 2.** We now show that for all $x \in [0, s_n)$, there holds

$$xg'_{n}(x)/n + (2 - 1/n)g_{n}(x) = g_{n}^{2}(x) + 2x. \qquad (3.10)$$
Write $g_n(x) = h_n(2n\sqrt{2x})/(2n)$, where $h_n(y) = yJ_{2n}(y)/J_{2n-1}(y)$. Observing that $h_n(y) = (y^{2n}J_{2n}(y))/(y^{2n-1}J_{2n-1}(y))$ and using (3.7) and then (3.6),

$$ 
\begin{aligned}
h'_{n}(y) &= \frac{y^{2n} J_{2n-1}(y)}{y^{2n-1} J_{2n-1}(y)} - \frac{y^{2n} J_{2n}(y) y^{2n-1} J_{2n-2}(y)}{(y^{2n-1} J_{2n-1}(y))^2} \\
&= y - h_n(y) \frac{J_{2n-2}(y)}{J_{2n-1}(y)} \\
&= y - h_n(y) \left[ \frac{-J_{2n}(y)}{J_{2n-1}(y)} + \frac{2(2n-1)}{y} \right] \\
&= y + \frac{h_n^2(y)}{y} - 2(2n-1) \frac{h_n(y)}{y}.
\end{aligned}
$$

But $g'_n(x) = h'_n(2n\sqrt{2x})/\sqrt{2x}$, whence

$$ 
\begin{aligned}
\frac{xg'_{n}(x)}{n} &= \frac{x}{n\sqrt{2x}} \left[ 2n\sqrt{2x} + \frac{h_{n}^{2}(2n\sqrt{2x})}{2n\sqrt{2x}} - 2(2n-1)\frac{h_{n}(2n\sqrt{2x})}{2n\sqrt{2x}} \right] \\
&= 2x + g_{n}^{2}(x) - (2 - 1/n)g_{n}(x).
\end{aligned}
$$

Step 3. Let us check that the sequence $(\beta_k^n)_{k \ge 1}$ satisfies (3.1). We already know that $\beta_1^n = 1$. Using (3.10),

$$ \sum_{k \ge 1} k(\beta_k^n/n)x^k + \sum_{k \ge 1} (2 - 1/n)\beta_k^n x^k = \sum_{k \ge 2} x^k \sum_{i=1}^{k-1} \beta_i^n \beta_{k-i}^n + 2x. $$

Thus for all $k \ge 2$, $(k/n)\beta_k^n + (2-1/n)\beta_k^n = \sum_{i=1}^{k-1} \beta_i^n \beta_{k-i}^n$ as desired. Consequently, $(\beta_k^n)_{k \ge 1} = (\alpha_k^n)_{k \ge 1}$, whence $f_n = g_n$ and $r_n = s_n = (j_{2n-1}/(2n))^2/2$.

Step 4. We know that $f_n$ is $C^\infty$ and increasing on $[0, r_n)$, that $f_n(0) = 0$ and that $\lim_{x \to r_n} f_n(x) = \infty$ (due to (3.9) and because $j_{2n} > j_{2n-1}$). Hence, there exists a unique $q_n \in [0, r_n)$ such that $f_n(q_n) = 1$ and we have $0 < f'_n(q_n) < \infty$. Applying (3.10) at $x = q_n$ (recall that $f_n = g_n$), we deduce that $q_n f'_n(q_n)/n + (2 - 1/n) = 1 + 2q_n$, so that $f'_n(q_n) = (n/q_n)[2q_n - 1 + 1/n]$.

Step 5. Put $\gamma_n = 2n\sqrt{2q_n} \in (0, 2n\sqrt{2r_n}) = (0, j_{2n-1})$. Then $f_n(q_n) = 1$ may be rewritten $\gamma_n J_{2n}(\gamma_n) = 2n J_{2n-1}(\gamma_n)$. We now prove that $j'_{2n-1} \le \gamma_n \le j'_{2n}$.

• First, using (3.5) with $k = 2n - 1$, we get $xJ_{2n}(x) = (2n - 1)J_{2n-1}(x) - xJ'_{2n-1}(x)$. Since $\gamma_n J_{2n}(\gamma_n) = 2nJ_{2n-1}(\gamma_n)$, we find $J_{2n-1}(\gamma_n) + \gamma_n J'_{2n-1}(\gamma_n) = 0$. Thus $\gamma_n \ge j'_{2n-1}$, because for $0 < x < j'_{2n-1} < j_{2n-1}$, $J_{2n-1}(x)$ and $J'_{2n-1}(x)$ are positive.

• We next show that $\gamma_n \le j'_{2n}$. We already know that $\gamma_n < j_{2n-1}$. We thus assume below that $j'_{2n} < j_{2n-1}$, because else, there is nothing to do. There holds $(2n/\gamma_n)[\gamma_n^2 J_{2n-1}(\gamma_n)]/[\gamma_n^2 J_{2n}(\gamma_n)] = 1$, whence $(2n/\gamma_n)[\log(\gamma_n^2 J_{2n}(\gamma_n))]' = 1$ by (3.7). Consequently,

$$ [\log J_{2n}(\gamma_n)]' = (\gamma_n/2n) - (2n/\gamma_n). $$

We know by (3.8) that $j'_{2n} > 2n$. Hence $J'_{2n} > 0$ and thus $[\log J_{2n}]' > 0$ on $[0, 2n]$. Thus $\gamma_n > 2n$, because for $x \le 2n$, we have $[\log J_{2n}(x)]' > 0$ and $(x/2n) - (2n/x) \le 0$. Hence $(\gamma_n/2n) - (2n/\gamma_n) > 0$. Thus $\gamma_n < j'_{2n}$, because for $x \in (j'_{2n}, j_{2n-1})$, $J'_{2n}(x) < 0$ so that $[\log J_{2n}(x)]' < 0$. Here we used that the second zero $j'_{2n,2}$ of
$J'_{2n}$ is greater than $j_{2n-1}$. Indeed, we have $j'_{2n,2} \ge j_{2n}$ (see Finch (2003)) and $j_{2n} \ge j_{2n-1}$ as already mentioned.

We thus have checked that $\gamma_n \in (j'_{2n-1}, j'_{2n})$. Using (3.8), we deduce that

$$ (2n - 1) + 2^{-1/3}|a'_{1}|(2n - 1)^{1/3} < \gamma_n < 2n + 2^{-1/3}|a'_{1}|(2n)^{1/3} + C(2n)^{-1/3}, $$

so that $\gamma_n = 2n + |a'_{1}|n^{1/3} + O(1)$. Recalling that $q_n = (\gamma_n/(2n))^2/2$, we easily deduce that $q_n = (1 + |a'_{1}|n^{-2/3})/2 + O(n^{-1})$ as desired. $\square$

We now have all the tools to give the

*Proof of Theorem 2.2-(i).* Fix $n \in \mathbb{N}$. Due to Lemma 3.2, there is a unique $q_n > 0$ such that $f_n(q_n) = 1$ and $f'_n(q_n) < \infty$. Applying Lemma 3.1, we deduce the existence and uniqueness of a solution to $(E_n)$. $\square$

Let us now give two weak forms of the equilibrium equation $(E_n)$.

**Lemma 3.3.** For $n \in \mathbb{N}$, consider the unique solution $(c_k^n)_{k \ge 1}$ to $(E_n)$. For any $\phi, \psi: \mathbb{N} \mapsto \mathbb{R}$ with at most polynomial growth,

$$ \frac{1}{m_0(c^n)} \sum_{k,l \ge 1} [\phi(k+l) - \phi(k) - \phi(l)] c_k^n c_l^n = \frac{1}{n} \sum_{k \ge 2} [\phi(k) - k\phi(1)] (k-1)c_k^n, \quad (3.11) $$

$$ \frac{2}{m_0(c^n)} \sum_{k,l \ge 1} [\psi(k+l) - \psi(k)] k c_k^n c_l^n = \frac{1}{n} \sum_{k \ge 2} [\psi(k) - \psi(1)] k(k-1)c_k^n. \quad (3.12) $$

*Proof:* Let $n \in \mathbb{N}$ be fixed.

**Step 1.** We first check that there is $u_n > 1$ such that $\sum_{k \ge 1} u_n^k c_k^n < \infty$. This allows us to justify the convergence of all the series in Steps 2 and 3 below. We know (see Lemma 3.1) that $c_k^n = \alpha_k^n q_n^{k-1}/f'_n(q_n)$. Hence for any $u > 0$, $\sum_{k \ge 1} u^k c_k^n = f_n(q_n u)/(q_n f'_n(q_n))$. Since the radius of convergence $r_n$ of the power series $f_n$ satisfies $r_n > q_n$ by Lemma 3.2, the result follows: choose $u_n > 1$ such that $q_n u_n < r_n$.

**Step 2.** We now prove (3.11). Multiply (2.3)-(2.2) by $\phi(k)$ and sum for $k \ge 1$. We get

$$ \frac{1}{m_0(c^n)} \left[ \sum_{k \ge 2} \phi(k) \sum_{i=1}^{k-1} c_i^n c_{k-i}^n - 2m_0(c^n) \sum_{k \ge 1} \phi(k) c_k^n \right] = \frac{1}{n} \left[ \sum_{k \ge 2} \phi(k)(k-1)c_k^n - [m_2(c^n) - 1]\phi(1) \right]. $$

But $\sum_{k \ge 2} \phi(k) \sum_{i=1}^{k-1} c_i^n c_{k-i}^n = \sum_{i,j \ge 1} \phi(i+j)c_i^n c_j^n$. Furthermore, there holds $2m_0(c^n) \sum_{k \ge 1} \phi(k)c_k^n = \sum_{i,j \ge 1} (\phi(i)+\phi(j))c_i^n c_j^n$, as well as $[m_2(c^n) - 1]\phi(1) = [m_2(c^n) - m_1(c^n)]\phi(1) = \sum_{k \ge 1} k\phi(1)(k-1)c_k^n$. This ends the proof of (3.11).

**Step 3.** To check (3.12), it suffices to apply (3.11) to the function $\phi(k) = k\psi(k)$ and to use that by symmetry, $\sum_{k,l \ge 1} [(k+l)\psi(k+l) - k\psi(k) - l\psi(l)]c_k^n c_l^n = 2 \sum_{k,l \ge 1} k[\psi(k+l) - \psi(k)]c_k^n c_l^n$. $\square$

**Remark 3.4.** For $n \in \mathbb{N}$, consider the solution $(c_k^n)_{k \ge 1}$ to $(E_n)$, write $p_k^n = c_k^n/m_0(c^n)$ and then $d_k^n = p_k^n/k$. Then one easily checks, starting from (3.11), that the sequence
$(d_k^n)_{k \ge 1}$ solves, for all $\phi: \mathbb{N} \to \mathbb{R}$ with at most polynomial growth,

$$ \frac{1}{2} \sum_{k,l \ge 1} [\phi(k+l) - \phi(k) - \phi(l)] k l d_k^n d_l^n = \sum_{k \ge 1} [\phi(k) - k\phi(1)] \frac{k(k-1)}{2n} d_k^n $$

and has total mass $\sum_{k \ge 1} k d_k^n = 1$. Hence $(d_k^n)_{k \ge 1}$ is an equilibrium state for a coagulation-fragmentation equation with multiplicative coagulation kernel (particles with masses $k, l$ merge at rate $kl$) and where each particle with mass $k$ breaks into $k$ particles with mass 1 at rate $k(k-1)/(2n)$. This explains the similarities between Theorem 2.3 and the results found by Ráth and Tóth (2009): the (non size-biased) particles' mass distribution $(p_k^n)_{k \ge 1}$ has the same limit, as $n \to \infty$, as the size-biased particles' mass distribution of the model considered by Ráth-Tóth, see Ráth and Tóth (2009, Eq. (15)). Observe however that the fragmentation rate in Ráth and Tóth (2009) is rather $(k-1)/n$, which thus differs from $k(k-1)/(2n)$. What seems important is just that roughly, for $n$ very large, only huge particles break down into atoms.

We are ready to handle the

*Proof of Theorem 2.2-(ii).* First, we know from Lemmas 3.1 and 3.2 that

$$ m_0(c^n) = \frac{1}{q_n f'_n(q_n)} = \frac{1}{n(2q_n - 1 + 1/n)} = \frac{1}{n[|a'_1|n^{-2/3} + O(1/n)]} \sim \frac{1}{|a'_1|n^{1/3}} $$

Next we use (3.11) with $\phi(k) = -1$. This gives

$$ m_0(c^n) = \frac{m_2(c^n) + m_0(c^n) - 2m_1(c^n)}{n}, \qquad (3.13) $$

whence, since $m_1(c^n) = 1$,

$$ m_2(c^n) = (n-1)m_0(c^n) + 2 \sim \frac{n^{2/3}}{|a'_1|} $$

Theorem 2.2-(ii) is established. $\square$

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

We can now prove the convergence of the particles' mass distribution.

*Proof of Theorem 2.3.* Let us put, for $k \ge 1$,

$$ p_k = \frac{2}{4^k k} \binom{2k-2}{k-1} $$

Using the Stirling formula, one immediately checks (2.5). Next, recall that the Catalan numbers, defined by $C_i = \frac{1}{i+1} \binom{2i}{i}$ for all $i \ge 0$, satisfy, see Richard and Weisstein (2014)

$$ C_i = \sum_{j=0}^{i-1} C_j C_{i-1-j}, \quad (i \ge 1), \\ \sum_{i \ge 0} C_i x^i = (1 - \sqrt{1 - 4x})/(2x), \quad x \in [0, 1/4]. $$

Observing that $p_k = 2 \cdot 4^{-k} C_{k-1}$, we easily deduce that $p_k = \frac{1}{2} \sum_{i=1}^{k-1} p_i p_{k-i}$ for $k \ge 2$, as well as $\sum_{k \ge 1} p_k = 2 \sum_{k \ge 1} C_{k-1}(1/4)^k = [\sum_{i \ge 0} C_i(1/4)^i]/2 = 1$.

We now check that (2.4) has at most one nonnegative solution. To this end, it suffices to show that (2.4) implies that $p_1 = 1/2$ (because this will determine the
value of $p_2 = (1/2)p_1^2$, of $p_3 = (1/2)[p_1p_2 + p_2p_1]$ and so on). To this end, it suffices
to write $p_1 = 1 - \sum_{k \ge 2} p_k = 1 - (1/2) \sum_{k \ge 2} \sum_{i=1}^{k-1} p_i p_{k-i} = 1 - (1/2)(\sum_{k \ge 1} p_k)^2 =
1 - 1/2 = 1/2$.

It only remains to prove that $\lim_n \sum_{k \ge 1} |p_k^n - p_k| = 0$. Since $\sum_{k \ge 1} p_k^n = \sum_{k \ge 1} p_k = 1$, it classically suffices to prove that $\lim_n p_k^n = p_k$ for all $k \ge 1$. First of all, we observe from (2.3) and Theorem 2.2-(ii) that

$$
p_1^n = \frac{c_1^n}{m_0(c^n)} = \frac{m_2(c^n) - 1}{2nm_0(c^n)} \sim \frac{n^{2/3}/|a'_1|}{2n/(|a'_1|n^{1/3})} \to \frac{1}{2} = p_1
$$

as $n \to \infty$. Next, we work by induction on $k$. Assume thus that for some $k \ge 2$,
$\lim_n p_l^n = p_l$ for $l = 1, \dots, k-1$. Then, using (2.2), we deduce that

$$
p_k^n = \frac{c_k^n}{m_0(c^n)} = \frac{1}{2 + (k-1)/n} \sum_{i=1}^{k-1} p_i^n p_{k-i}^n \to \frac{1}{2} \sum_{i=1}^{k-1} p_i p_{k-i} = p_k
$$

as $n \to \infty$. This concludes the proof.

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

We now study the size-biased particles’ mass distribution. We start with the
computation of all the moments and deduce a convergence result.

**Lemma 3.5.** For each $n \in \mathbb{N}$, consider the unique solution $(c_k^n)_{k \ge 1}$ to $(E_n)$ and
the probability measure $\mu_n = \sum_{k \ge 1} k c_k^n \delta_{n^{-2/3}k}$ on $(0, \infty)$.

(i) For any $i \ge 1$,

$$
m_{i+1}(c^n) \xrightarrow{n \to \infty} M_i n^{2i/3},
$$

where the sequence $(M_i)_{i \ge 0}$ is defined by $M_0 = 1$, $M_1 = 1/|a'_1|$ and, for $i \ge 1$,

$$
M_{i+1} = 2|a'_1| \sum_{j=0}^{i-1} \binom{i}{j} M_j M_{i-j-1}.
$$

(ii) There is a probability measure $\mu$ on $[0, \infty)$ such that for any $\phi \in C([0, \infty))$
with at most polynomial growth,

$$
\sum_{k \ge 1} \phi(n^{-2/3}k) k c_k^n = \int_0^\infty \phi(x) \mu_n(dx) \xrightarrow{n \to \infty} \int_0^\infty \phi(x) \mu(dx).
$$

This probability measure satisfies, for all $i \ge 0$, $\int_0^\infty x^i \mu(dx) = M_i$.

Let us mention that the recursive formula for the moment sequence $(M_i)_{i \ge 0}$
resembles that due to Flajolet and Louchard for the moments of the area of the
Brownian excursion, see Equation (9) in Janson (2007).

*Proof:* Applying (3.12) with $\psi(k) = k^i$ with $i \ge 1$, we easily get

$$
\frac{2}{m_0(c^n)} \sum_{j=0}^{i-1} \binom{i}{j} m_{j+1}(c^n) m_{i-j}(c^n) = \frac{1}{n} (m_{i+2}(c^n) - m_{i+1}(c^n) - m_2(c^n) + 1),
$$

so that

$$
m_{i+2}(c^n) = m_{i+1}(c^n) + m_2(c^n) - 1 + \frac{2n}{m_0(c^n)} \sum_{j=0}^{i-1} \binom{i}{j} m_{j+1}(c^n) m_{i-j}(c^n).
$$

From this and the fact that we already know that $m_2(c^n) \sim M_1 n^{2/3}$ and that
$m_0(c^n) \sim n^{-1/3}/|a'_1|$, one can easily check point (i) by induction.
Next, it holds that $\int_0^\infty x^i \mu_n(dx) = n^{-2i/3} m_{i+1}(c^n)$ for all $i \ge 1$. Consequently, we know from point (i) that $\lim_{n \to \infty} \int_0^\infty x^i \mu_n(dx) = M_i$ for any $i \ge 1$.

But there is a unique probability measure $\mu$ on $[0, \infty)$ such that $\int_0^\infty x^i \mu(dx) = M_i$ for all $i \ge 1$. Indeed, an immediate (and rough) induction using that $|a'_1| > 1$ shows that $M_i \le (2|a'_1|)^{i!}$ for all $i \ge 1$. This implies the finiteness of an exponential moment for $\mu$, so that $\mu$ is characterized by its moments. As a conclusion, $\mu_n$ goes weakly, as $n \to \infty$, to $\mu$. This shows point (ii) for all $\phi: [0, \infty) \mapsto \mathbb{R}$ continuous and bounded. The extension to continuous functions with at most polynomial growth easily follows from point (i). $\square$

Let us now show that $\mu$ satisfies some equilibrium equation.

**Lemma 3.6.** Consider the probability measure $\mu$ on $[0, \infty)$ defined in Lemma 3.5. For all $\phi \in C^1([0, \infty))$ such that $\phi$ and $\phi'$ have at most polynomial growth,

$$2|a'_1| \int_0^\infty \int_0^\infty \left[ \frac{\phi(x+y) - \phi(x)}{y} \mathbb{1}_{\{y>0\}} + \phi'(x) \mathbb{1}_{\{y=0\}} \right] \mu(dx)\mu(dy) \\ = \int_0^\infty x[\phi(x) - \phi(0)]\mu(dx).$$

*Proof:* We consider $\phi$ as in the statement and, for $n \in \mathbb{N}$, $\mu_n = \sum_{k \ge 1} k c_k^n \delta_{n^{-2/3}k}$ as in Lemma 3.5. Apply (3.12) with $\psi(k) = \phi(n^{-2/3}k)$:

$$\frac{2n^{-2/3}}{m_0(c^n)} \sum_{k,l \ge 1} \frac{\phi(n^{-2/3}(k+l)) - \phi(n^{-2/3}k)}{n^{-2/3}l} k c_k^n l c_l^n \\ = \frac{n^{2/3}}{n} \sum_{k \ge 1} (\phi(n^{-2/3}k) - \phi(n^{-2/3})) (n^{-2/3}k - n^{-2/3}) k c_k^n.$$

Multiply this equality by $n^{1/3}$. In terms of $\mu_n$, this can be written as

$$\frac{2n^{-1/3}}{m_0(c^n)} \int_0^\infty \int_0^\infty \frac{\phi(x+y) - \phi(x)}{y} \mu_n(dx)\mu_n(dy) \\ = \int_0^\infty [\phi(x) - \phi(n^{-2/3})](x - n^{-2/3})\mu_n(dx).$$

Recall that $m_0(c^n) \sim n^{-1/3}/|a'_1|$ as $n \to \infty$, so that

$$2|a'_1| \int_0^\infty \int_0^\infty \frac{\phi(x+y) - \phi(x)}{y} \mu_n(dx)\mu_n(dy) \\ \sim \int_0^\infty [\phi(x) - \phi(n^{-2/3})](x - n^{-2/3})\mu_n(dx) \\ \sim \int_0^\infty [\phi(x) - \phi(0)]x\mu_n(dx).$$

To obtain the last equivalent, use that $\phi$ is continuous, that $\sup_n \int_0^\infty x \mu_n(dx) < \infty$ and that $\sup_n \int_0^\infty |\phi(x) - \phi(0)| \mu_n(dx) < \infty$ by Lemma 3.5-(i) since $\phi$ has at most polynomial growth. Define now the function $\Gamma(x, y) = ([\phi(x+y) - \phi(x)]/y) \mathbb{1}_{\{y>0\}} + \phi'(x) \mathbb{1}_{\{y=0\}}$, which is continuous and has at most polynomial growth on $[0, \infty)^2$.
Since $\mu_n$ does not give weight to 0, we have

$$ 
\begin{aligned}
& 2|a'_1| \int_0^\infty \int_0^\infty \Gamma(x, y) \mu_n(dx) \mu_n(dy) && (3.14) \\
&= 2|a'_1| \int_0^\infty \int_0^\infty \frac{\phi(x+y) - \phi(x)}{y} \mu_n(dx) \mu_n(dy) \\
&\sim \int_0^\infty [\phi(x) - \phi(0)] x \mu_n(dx).
\end{aligned}
$$

Recall that for any $\psi \in C([0, \infty))$ with at most polynomial growth, there holds $\int_0^\infty \psi(x)\mu_n(dx) \xrightarrow{n \to \infty} \int_0^\infty \psi(x)\mu(dx)$ due to Lemma 3.5. This implies that for all $\Psi \in C([0, \infty)^2)$ with at most polynomial growth, $\int_0^\infty \int_0^\infty \Psi(x, y)\mu_n(dx)\mu_n(dy) \xrightarrow{n \to \infty} \int_0^\infty \int_0^\infty \Psi(x, y)\mu(dx)\mu(dy)$. Taking the limit as $n \to \infty$ in (3.14), we deduce that

$$ 2|a'_1| \int_0^\infty \int_0^\infty \Gamma(x, y)\mu(dx)\mu(dy) = \int_0^\infty [\phi(x) - \phi(0)]x\mu(dx) $$

as desired.

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

We now try to determine a quantity resembling the Laplace transform of $\mu$. This is a usual trick for Smoluchowski's coagulation with constant kernel, see e.g. Deaconu and Tanré (2000).

**Lemma 3.7.** Consider the probability measure $\mu$ defined in Lemma 3.5. Then for all $q \ge 0$,

$$ \ell(q) := \int_0^\infty (1 - e^{-qx}) \frac{\mu(dx)}{x} = \frac{-\mathrm{Ai}'(q + a'_1)}{|a'_1|\mathrm{Ai}(q + a'_1)}. $$

*Proof:* We put $\beta = \mu(\{0\})$, which we cannot exclude to be nonzero at the moment. We apply Lemma 3.6 with $\phi(x) = (1 - e^{-qx})/x$, which is indeed in $C^1([0, \infty))$ with $\phi(0) = q$ and $\phi'(0) = -q^2/2$. This yields

$$ 2|a'_1|[A_1(q) + A_2(q) + A_3(q) + A_4(q)] = B_1(q) - B_2(q), $$

where

$$ 
\begin{aligned}
A_1(q) &= \int_0^\infty \int_0^\infty \frac{\phi(x+y) - \phi(x)}{y} \mathbf{1}_{\{x>0, y>0\}} \mu(dx)\mu(dy), \\
A_2(q) &= \beta \int_0^\infty \frac{\phi(y) - \phi(0)}{y} \mathbf{1}_{\{y>0\}} \mu(dy), \\
A_3(q) &= \beta \int_0^\infty \phi'(x) \mathbf{1}_{\{x>0\}} \mu(dx), \\
A_4(q) &= \beta^2 \phi'(0), \\
B_1(q) &= \int_0^\infty x\phi(x)\mu(dx), \\
B_2(q) &= \phi(0) \int_0^\infty x\mu(dx).
\end{aligned}
$$

Recalling Lemma 3.5 and that $\phi(0) = q$, we see that $B_2(q) = q/|a'_1|$. Next,

$$ B_1(q) = \int_0^\infty (1 - e^{-qx}) \mu(dx) = 1 - \ell'(q). $$
One can check that $\phi'(0) = -q^2/2$, whence $A_4(q) = -\beta^2 q^2/2$. A computation shows that $[\phi(x) - \phi(0)]/x + \phi'(x) = -q\phi(x)$, so that

$$A_2(q) + A_3(q) = -\beta q \int_0^\infty x^{-1}(1 - e^{-qx}) \mathbb{1}_{\{x>0\}} \mu(dx) = -\beta q[\ell(q) - \beta q].$$

Finally, using a symmetry argument and then that $(x+y)\phi(x+y) - x\phi(x) - y\phi(y) = -xy\phi(x)\phi(y)$,

$$\begin{aligned}
A_1(q) &= \int_0^\infty \int_0^\infty \frac{x\phi(x+y) - x\phi(x)}{xy} \mathbb{1}_{\{x>0, y>0\}} \mu(dx)\mu(dy) \\
&= \frac{1}{2} \int_0^\infty \int_0^\infty \frac{(x+y)\phi(x+y) - x\phi(x) - y\phi(y)}{xy} \mathbb{1}_{\{x>0, y>0\}} \mu(dx)\mu(dy) \\
&= -\frac{1}{2} \left( \int_0^\infty (1 - e^{-qx}) \mathbb{1}_{\{x>0\}} \frac{\mu(dx)}{x} \right)^2 \\
&= -\frac{1}{2} (\ell(q) - \beta q)^2.
\end{aligned}$$

All this shows that

$$2|a'_1| \left[ -\frac{1}{2}(\ell(q) - \beta q)^2 - \beta q[\ell(q) - \beta q] - \beta^2 q^2/2 \right] = 1 - \ell'(q) - q/|a'_1|,$$

whence

$$\ell'(q) = 1 - q/|a'_1| + |a'_1|\ell^2(q). \qquad (3.15)$$

This equation, together with the initial condition $\ell(0) = 0$, has a unique maximal solution due to the Cauchy-Lipschitz Theorem. And one can check that this unique maximal solution is nothing but

$$\ell(q) = \frac{-\mathrm{Ai}'(q + a'_1)}{|a'_1|\mathrm{Ai}(q + a'_1)},$$

which is defined for $q \in (a_1 - a'_1, \infty)$, because the Airy function does not vanish on $(a_1, \infty)$. Indeed, it suffices to use that since $\mathrm{Ai}''(x) = x\mathrm{Ai}(x)$,

$$\frac{d}{dx} \left( \frac{\mathrm{Ai}'(x)}{\mathrm{Ai}(x)} \right) = x - \left( \frac{\mathrm{Ai}'(x)}{\mathrm{Ai}(x)} \right)^2$$

and that by definition, $\mathrm{Ai}'(a'_1) = 0$. $\square$

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

We now write down two formulae of Darling (1983) and Louchard (1984) that we found in the survey paper of Janson (2007, p 94). Denote by $(e_t)_{t \in [0,1]}$ the normalized Brownian excursion and define its area as $B_{\mathrm{ex}} = \int_0^1 e_t dt$. Put, for $y \ge 0$,

$$\psi_{\mathrm{ex}}(y) = \mathbb{E}[e^{-yB_{\mathrm{ex}}}] \qquad (3.16)$$

There hold

$$\psi_{\mathrm{ex}}(y) = \sqrt{2\pi y} \sum_{j=1}^{\infty} \exp(-2^{-1/3}|a_j|y^{2/3}), \qquad (3.17)$$

$$\int_0^\infty (1 - e^{-qy}) \frac{\psi_{\mathrm{ex}}(y^{3/2})}{\sqrt{2\pi y^3}} dy = 2^{1/3} \left( \frac{\mathrm{Ai}'(0)}{\mathrm{Ai}(0)} - \frac{\mathrm{Ai}'(2^{1/3}q)}{\mathrm{Ai}(2^{1/3}q)} \right). \qquad (3.18)$$

This allows us to find a link between our probability measure $\mu$ and $\psi_{\mathrm{ex}}$.
**Lemma 3.8.** Consider the probability measure $\mu$ defined in Lemma 3.5. Then

$$ \frac{\mu(dx)}{x} = \frac{\psi_{\text{ex}}(\sqrt{2}x^{3/2})e^{|a'_1|x}}{2\sqrt{\pi}|a'_1|x^{3/2}} \mathbf{1}_{\{x>0\}}dx. $$

*Proof:* Using Lemma 3.7, we deduce that for all $q \ge 0$,

$$ \begin{aligned} \int_0^\infty (1 - e^{-qx}) \frac{e^{-|a'_1|x} \mu(dx)}{x} &= \int_0^\infty (e^{-|a'_1|x} - 1) \frac{\mu(dx)}{x} + \int_0^\infty (1 - e^{-(|a'_1|+q)x}) \frac{\mu(dx)}{x} \\ &= -\ell(|a'_1|) + \ell(q + |a'_1|) \\ &= \frac{1}{|a'_1|} \left( \frac{\text{Ai}'(0)}{\text{Ai}(0)} - \frac{\text{Ai}'(q)}{\text{Ai}(q)} \right). \end{aligned} $$

Next, (3.18) implies that for all $q \ge 0$,

$$ \begin{aligned} \int_0^\infty (1 - e^{-qx}) \frac{\psi_{\text{ex}}(\sqrt{2}x^{3/2})}{2\sqrt{\pi}|a'_1|x^{3/2}} dx &= \int_0^\infty (1 - e^{-q2^{-1/3}y}) \frac{\psi_{\text{ex}}(y^{3/2})}{|a'_1|\sqrt{2\pi y^3}} 2^{-1/3} dy \quad (3.19) \\ &= \frac{1}{|a'_1|} \left( \frac{\text{Ai}'(0)}{\text{Ai}(0)} - \frac{\text{Ai}'(q)}{\text{Ai}(q)} \right). \end{aligned} $$

We conclude by injectivity of the Laplace transform. □

We may finally give the

*Proof of Theorem 2.4.* Define $c: (0, \infty) \mapsto (0, \infty)$ by

$$ c(x) = \frac{\psi_{\text{ex}}(\sqrt{2}x^{3/2})e^{|a'_1|x}}{2\sqrt{\pi}|a'_1|x^{3/2}}. \quad (3.20) $$

By Lemma 3.8, the probability measure $\mu$ defined in Lemma 3.5 is nothing but $\mu(dx) = xc(x)dx$. We thus have $\int_0^\infty xc(x)dx = 1$. Recalling Lemma 3.5, we know that $\lim_{n \sum_{k \ge 1} \phi(n^{-2/3}k)} kc_n^k = \int_0^\infty \phi(x)xc(x)dx$ for any $\phi \in C([0, \infty))$ with at most polynomial growth. Using (3.17), we immediately deduce that

$$ c(x) = |a'_1|^{-1} e^{|a'_1| x} \sum_{j=1}^{\infty} e^{-|a_j| x}. \quad (3.21) $$

It is clear from (3.21) that $c \in C^\infty((0, \infty))$ and that $c(x) \stackrel{\sim}{\sim} |a'_1|^{-1}e^{|a'_1|-|a_1|}x$: it suffices to use that $0 < |a_1| < |a_2| < \dots$, that $|a_j| \stackrel{j \to \infty}{\sim} (3\pi j/2)^{2/3}$ and Lebesgue's dominated convergence theorem. It is immediate from (3.20) that $c(x) \stackrel{\sim}{\sim} 1/(2\sqrt{\pi}|a'_1|x^{3/2})$. Since we now know that $\mu(dx) = xc(x)dx$ does not give weight to zero, (2.6) follows from Lemma 3.6. Finally, (2.7) follows from Lemma 3.7 when $q \ge 0$. It is easily extended to $q > a_1 - a'_1$ using that $c(x) \stackrel{\sim}{\sim} |a'_1|^{-1}e^{|a'_1|-|a_1|}x$. □

# 4. Perfect simulation algorithms

In this section, we provide some perfect simulation algorithms: we introduce a pruning procedure $\mathcal{P}_n$ of trees (for each $n \ge 1$) that will allow us to interpret

* the particles' mass distribution $(p_k^n)_{k \ge 1}$ as the law of the number of leaves of $G_n = \mathcal{P}_n(G)$, where $G$ is a binary critical Galton-Watson tree;
• the size-biased particles' mass distribution $(kc_k^n)_{k \ge 1}$ as the law of the number of leaves of $\hat{G}_n = \mathcal{P}_n(\hat{G})$, where $\hat{G}$ is the so-called size-biased binary critical Galton-Watson tree defined in Lyons et al. (1995, Section 5).

We also show that $G_n$ obviously tends to $G$ as $n \to \infty$, which gives a probabilistic interpretation of Theorem 2.3. Let us insist on the fact that the goal of this section is definitely not numerical: we only want to provide some probabilistic interpretations.

**4.1. The pruning procedure.** Consider a rooted discrete binary tree $T$, that is a set of vertices $V(T)$ and of edges $E(T)$, satisfying the usual properties of binary trees. The root is denoted by $\emptyset$ and we add a vertex $\star$ and an edge joining $\star$ to $\emptyset$. We denote by $L(T)$ the set of the leaves of $T$. We take the convention that $\star$ is not a leaf, but $\emptyset$ may be a leaf (iff $V(T)$ is reduced to $\{\star, \emptyset\}$). We denote by $N(T)$ the set of internal vertices (nodes) of $T$, that is $N(T) = V(T) \setminus (\{\star\} \cup L(T))$.

We recall that for any binary tree $T$, it holds that $|L(T)| = |N(T)| + 1$.

**Definition 4.1.** Let $T$ be a rooted discrete binary tree with at most one infinite branch and let $n \ge 1$. We define the (random) subtree $\mathcal{P}_n(T)$ as follows.

**Step 1.** We now think of each edge $e$ as a line segment with the length $\kappa_e$, where $(\kappa_e)_{e \in E(T)}$ is an i.i.d. family of $\mathcal{Exp}(2)$-distributed random variables. This induces a distance $d$ on $T$. For each $t \ge 0$, define $T_t = \{x \in T : d(x, \star) = t\}$ and let $t_0 = \inf\{t > 0 : T_t = \emptyset\} \le \infty$ be the height of $T$.

**Step 2.** For each internal node $z \in N(T)$, consider the branch $B_T(z)$ (endowed with the lengths introduced at Step 1) joining $z$ to $\star$ and consider a Poisson point process $\pi_z^n$ with rate $1/n$ on $B_T(z)$ (conditionally on $(\kappa_e)_{e \in E(T)}$, these Poisson processes are taken mutually independent). All the marks of these Poisson process are activated. In words, we will say that the marks of $\pi_z^n$ are sent by the node $z$ on the branch $B_T(z)$.

**Step 3.** We explore the tree from the top until we arrive at $\star$ (that is, we consider $T_t$ for $t$ decreasing from $t_0$ to 0) with the following rule: each time we encounter an active mark of a Poisson point process (defined in Step 2), we remove all the subtree above the mark (and replace it by a leaf) and we deactivate all the marks sent by nodes in this subtree. See Figure 4.1 for an illustration.

**Step 4.** We call $\mathcal{P}_n(T)$ the resulting tree, which is a subtree of $T$ containing $\star$.

A mark of a Poisson process is said to be useful if it has generated a pruning at some time in the procedure and useless otherwise.

Let us insist on the fact that we explore the tree from the top (i.e., starting from points that are far away from the root), but this exploration order concerns the *cut-points* and not the *nodes* that send them. For example in Figure 4.1, the cut-point sent by $x$ will be explored before the one sent by $y$, even if $y$ is higher (further from the root).

If $T$ is infinite, then Step 3 looks ill-posed at first glance, since the top of the tree lies at infinity.

**Remark 4.2.** Let $T$ be a tree with one infinite branch. Then Definition 4.1 makes sense. Indeed, consider the first (starting from $\star$) edge $e = (x_1, x_2)$ on this infinite branch, with $x_2$ child of $x_1$, such that $\pi_{x_2}^n$ has a mark in $e$. This happens for each
edge independently with positive probability (not depending on $e$), so that such an edge a.s. exists. Then $e$ will be cut (either by $\pi_{x_2}^n$ or by a Poisson process corresponding to a node above $x_2$) and this will make inactive all the Poisson processes corresponding to nodes above $x_2$. Thus everything will happen as if $x_2$ was a leaf of $T$, neglecting all the marks due to Poisson processes corresponding to nodes above $x_2$. See Figure 4.2 for an illustration.

<figure><img src="image_6.png" /><figcaption>FIGURE 4.1. The pruning procedure $\mathcal{P}_n$</figcaption></figure>

On Figure A, the tree $T$ is drawn, its edges being endowed with i.i.d. $\mathcal{Exp}(2)$-distributed random variables. The marks of the Poisson processes are represented as follows: here $\pi_x^n$ has one mark (on the edge just under $x$), $\pi_y^n$ has one mark (on the edge above $\emptyset$ on the left), $\pi_z^n$ has three marks (one on the edge just under $z$, one on the edge above $\emptyset$ on the right, one on the edge $(\star, \emptyset)$) and $\pi_t^n$ has one mark on the edge just under $t$. All the other Poisson point processes have no mark.

Thus starting from the top of the tree, we first encounter the mark just under $z$. On Figure B, we have drawn the resulting tree: we have replaced the subtree above this mark by a leaf and we have erased (deactivated) the marks sent by nodes of this subtree (here, the two other marks sent by $z$).

Then we encounter the mark sent by $t$ and the resulting tree is drawn on Figure C. Finally, we encounter the mark sent by $x$, which makes inactive the mark sent by $y$ and the resulting tree is drawn on figure D. Since there are no marks any more, the tree of figure D is $\mathcal{P}_n(T)$.
4.2. Particles' mass distribution. We can now give an interpretation in terms of trees of the particles' mass distribution.

**Proposition 4.3.** Consider a binary critical Galton-Watson tree G (BCGWT in short), that is a Galton-Watson tree with offspring distribution $(\delta_0 + \delta_2)/2$. Fix $n \ge 1$ and let $G_n = \mathcal{P}_n(G)$, see Definition 4.1 (all the random objects used by $\mathcal{P}_n$ are taken conditionally on G).

(i) For all $k \ge 1$, $\Pr[|L(G)| = k] = p_k$, where $p_k$ was defined in Theorem 2.3.

(ii) For all $k \ge 1$, $\Pr[|L(G_n)| = k] = p_k^n$, where $p_k^n$ was defined in Theorem 2.3.

Clearly, $G_n = \mathcal{P}_n(G)$ can be perfectly simulated. We thus have a perfect simulation algorithm for $(p_k^n)_{k \ge 1}$.

*Heuristic proof.* Any particle with mass $k$ can be seen as a cluster of $k$ particles with mass 1. And it is natural to use a tree to represent the *genealogy* of this particle; if this particle has a mass $k$, then this tree will have $k$ leaves.

To be more precise, we need first to handle some computations. Divide (3.11) by $m_0(c^n)$, use that $\sum_{k,l \ge 1} \phi(l)p_k^n p_l^n = \sum_{k \ge 1} \phi(k)p_k^n$ and that, recalling (3.13),

$$ \frac{1}{n} \sum_{k \ge 2} (k-1)^2 \phi(1) p_k^n = \phi(1) \frac{m_2(c^n) + m_0(c^n) - 2m_1(c^n)}{n m_0(c^n)} = \phi(1) = \sum_{k \ge 1} \phi(1) p_k^n. $$

One gets, for all reasonable $\phi$,

$$ \sum_{k,l \ge 1} [\phi(k+l) - \phi(k)] p_k^n p_l^n + \sum_{k \ge 2} [\phi(1) - \phi(k)][1 + (k-1)/n] p_k^n = 0. $$

<figure><img src="image_7.png" /><figcaption>FIGURE 4.2. Applying $\mathcal{P}_n$ to a tree with one infinite branch</figcaption></figure>

The first edge $e = (x_1, x_2)$ on the infinite branch such that $\pi_{x_2}^n$ has a mark on $e$ is here the edge under $t$ and $x_2 = t$. Then independently of what can happen above $t$, this edge will be cut and all the marks sent by nodes above $t$ will be erased. Indeed, there are two possibilities, calling $M$ the mark sent by $t$.

Case 1: This mark $M$ is useful and thus $M$ is replaced by a leaf and all the marks sent between $\star$ and $M$ by nodes above $t$ are erased.

Case 2. This mark $M$ is useless and then it is necessarily deactivated by a useful mark lying between $M$ and $t$ (sent by a node above $t$). We conclude as in case 1.

In any case, it is not necessary to know what happens above $t$ to conclude that, with this configuration, $\mathcal{P}_n(T)$ will be the same as in Figure 4.1-D.
Thus $(p_k^n)_{k \ge 1}$ can be seen as the equilibrium of the mass of a particle with the following dynamics: (i) it merges with an independent similar particle at rate 1, (ii) its mass is reset to 1 at rate 1, (iii) its mass is reset to 1 at rate $(k-1)/n$, where $k$ is its mass.

Consider now such a particle at equilibrium. First neglect (iii) and follow the history of the particle backward in time: it has merged with a similar particle at rate 1 and it has been reset to 1 at rate 1. Thus this particle is subjected to events at rate 2 and each time an event occurs, it is a coalescence (node) with probability 1/2 and a breakage (leaf) with probability 1/2. This can be represented by a BCGWT, of which the edges have a length with law $\mathcal{E}xp(2)$. Now we take (iii) into account: we start from the past (i.e. from the top of the tree), we follow the branches of the tree and reset the particle to 1 at rate $(k-1)/n$, where $k$ is the mass of the particle, i.e. the number of leaves of the subtree above the point under consideration. Using finally that the number of nodes of a binary tree is precisely its number of leaves minus 1, we guess that $\mathcal{P}_n(G)$ should indeed provide a perfect simulation algorithm for the genealogy of a particle at equilibrium.

We now handle a rigorous proof.

*Proof:* We start with (i), which is completely standard. Define $q_k = \text{Pr}[|L(G)| = k]$ for $k \ge 1$ and use the branching property: knowing that $G$ is not reduced to the root, it can be written as $G = \{\emptyset\} \cup G' \cup G''$, for two independent copies $G'$, $G''$ of $G$. Hence conditionally on $\{|L(G)| \ge 2\}$, $L(G) = L(G') \cup L(G'')$, whence, for $k \ge 2$, $\{|L(G)| = k\} = \{|L(G)| \ge 2\} \cap \cup_{i=1}^{k-1}\{|L(G')| = i, |L(G'')| = k - i\}$ and thus $q_k = \frac{1}{2} \sum_{i=1}^{k-1} q_i q_{k-i}$. Next, it obviously holds that $q_1 = 1/2$. Since $p_1 = q_1$ and recalling (2.4), it follows that $q_k = p_k$ for all $k \ge 1$.

We now check (ii). Using the branching property of $G$ recalled above, it follows from the pruning procedure $\mathcal{P}_n$ that for $k \ge 2$,

$$ \{|L(G_n)| = k\} = \{|L(G)| \ge 2\} \cap \bigcup_{i=1}^{k-1} \{|L(G'_n)| = i, |L(G''_n)| = k - i, A_n\}, $$

where $A_n$ is the event that no pruning occurs in the edge $(\star, \emptyset)$ and where $G'_n$ and $G''_n$ are two independent copies of $G_n$, independent of $(\{|L(G)| \ge 2\}, (\pi_x^n|_{(\star, \emptyset)})_{x \in N(G)})$. Indeed, having a look at Figure 4.1, we see that if $|L(G)| \ge 2$, the pruning procedure, before exploring the edge $(\star, \emptyset)$, will produce two independent pruned Galton-Watson trees $G'_n$ and $G''_n$ (and it then only remains to explore the edge $(\star, \emptyset)$).

First observe that conditionally on $(G'_n, G''_n)$, the event $A_n$ occurs if for all $x \in N(G'_n) \cup N(G''_n) \cup \{\emptyset\}$, $\pi_x^n$ has no mark in $(\star, \emptyset)$. Recalling that the length of this edge is $\mathcal{E}xp(2)$-distributed, that the Poisson processes $\pi_x^n$ have rate $1/n$ and that here we have $|N(G'_n) \cup N(G''_n) \cup \{\emptyset\}| = |N(G'_n)| + |N(G''_n)| + 1 = |L(G'_n)| + |L(G''_n)| - 1$ Poisson processes, one easily deduces that

$$ \text{Pr}[A_n | G'_n, G''_n] = \frac{2}{(|L(G'_n)| + |L(G''_n)| - 1)/n + 2}. $$

Put now $r_k^n = \text{Pr}[|L(G_n)| = k]$. From the previous study, we get, for $k \ge 2$,

$$ r_k^n = \frac{1}{2} \sum_{i=1}^{k-1} r_i^n r_{k-i}^n \frac{2}{(i + k - i - 1)/n + 2}. $$
whence $[2+(k-1)/n]r_k^n = \sum_{i=1}^{k-1} r_i^n r_{k-i}^n$. Using the arguments and notation of the proof of Lemma 3.1-Step 1 (with $r=1$), we deduce that $r_k^n = \alpha_k^n (r_1^n)^k$ for all $k \ge 1$. But we also know, since $G_n \subset G$ is a.s. finite, that $\sum_{k \ge 1} \alpha_k^n (r_1^n)^k = \sum_{k \ge 1} r_k^n = 1$. Recalling (3.2) and Lemma 3.2, we conclude that $r_1^n = q_n$, whence $r_k^n = \alpha_k^n (q_n)^k$ for all $k \ge 1$. By Lemma 3.1, it also holds that $p_k^n = c_k^n/m_0(c^n) = \alpha_k^n(q_n)^k$, which concludes the proof. $\square$

**4.3. Probabilistic interpretation of Theorem 2.3.** This is not hard from Proposition 4.3. Clearly, as $n \to \infty$, the probability that $G = G_n$ tends to 1, implying that the law of $|L(G_n)|$, i.e. $(p_k^n)_{k \ge 1}$, tends to the law of $|L(G)|$, i.e. $(p_k)_{k \ge 1}$.

Indeed, we have $G_n = G$ as soon as for all $x \in N(G)$, $\pi_x^n$ has no mark on the branch $B_G(x)$ joining $\star$ to $x$ in $G$. Conditionally on $(G, (\kappa_e)_{e \in N(G)})$, this occurs with probability $\prod_{x \in N(G)} \exp(-\frac{1}{n} \sum_{e \in E(G), e \subset B_G(x)} \kappa_e)$, which a.s. tends to 1, since $G$ is a.s. a finite tree.

**4.4. Size-biased particles' mass distribution.** We now interpret the size-biased particles' mass distribution in terms of a pruned size-biased Galton-Watson tree.

**Definition 4.4.** Consider a family of i.i.d. binary critical Galton-Watson trees $(G(i))_{i \ge 1}$. We call size-biased binary critical Galton-Watson tree (SBBCGWT in short) the binary tree $\hat{G}$ with one infinite branch (called the backbone) as in Figure 4.3, where for each $i \ge 1$, we plant $G(i)$ on $\star i$.

<figure><img src="image_8.png" /><figcaption>FIGURE 4.3. A size-biased binary critical Galton-Watson tree</figcaption></figure>

Traditionally, each $G(i)$ is planted above or under the backbone at random, but this is absolutely useless for our purpose, since we never take into account any order on the vertices. See Lyons et al. (1995, Section 2) for some indications about the terminology *size-biased*.

**Proposition 4.5.** Consider a SBBCGWT $\hat{G}$ as in Definition 4.4 and fix $n \ge 1$. Let $\hat{G}_n = \mathcal{P}_n(\hat{G})$, recall Definition 4.1 (all the the random objects used by $\mathcal{P}_n$ are taken conditionally on $\hat{G}$). Then for all $k \ge 1$, $\Pr(|L(\hat{G}_n)| = k) = kc_k^n$, where $(c_k^n)_{k \ge 1}$ was defined in Theorem 2.2.

Since $\mathcal{P}_n(\hat{G})$ can be perfectly simulated due to Remark 4.2, this provides a perfect simulation algorithm for $(kc_k^n)_{k \ge 1}$. It seems striking that $\mathcal{P}_n(\hat{G})$ is a size-biased version of $\mathcal{P}_n(G)$, with $\hat{G}$ a size-biased version of $G$. However, the two notions of *size-biased* are quite different: $\mathcal{P}_n(\hat{G})$ is a version of $\mathcal{P}_n(G)$ *biased by the*
number of leaves, while $\hat{G}$ is a version of $G$ biased by the size of the population at generation $n$ (with $n \to \infty$).

*Heuristic proof.* First rewrite (3.12) as

$$ \sum_{k,l \ge 1} [\psi(k+l) - \psi(k)]kc_k^n 2p_l^n + \sum_{k \ge 2} [\psi(1) - \psi(k)] \frac{k-1}{n} kc_k^n = 0. $$

Thus $(kc_k^n)_{k \ge 1}$ can be seen as the equilibrium of the mass of a particle with the following dynamics: (i) it merges with an independent particle with law $(p_k^n)_{k \ge 1}$ at rate 2, (ii) its mass is reset to 1 at rate $(k-1)/n$, where $k$ is its mass.

Consider such a particle at equilibrium and first neglect (ii). Then obviously, we can represent its genealogy as a forest of (non size-biased) particles $P_n(G(i))$, the length of the edges one the backbone (those between these particles) being $\mathcal{E}xp(2)$-distributed. Then, take (ii) into account: start from the past and prune the edges on the backbone at rate $(k-1)/n$, where $k$ is the mass of the particle, i.e. the number of leaves of the subtree above the point under consideration, whence $k-1$ is the corresponding number of nodes.

What we really do is slightly different, since we prune the (non size-biased) particles and the backbone simultaneously, but one can easily get convinced that this changes nothing.

We now give some rigorous arguments.

*Proof:* Consider the problem with unknown $(t_k^n)_{k \ge 1}$ (here $(p_k^n)_{k \ge 1}$ is given)

$$ \sum_{k \ge 1} t_k^n = 1, \quad [(k-1)/n + 2]t_k^n = 2 \sum_{i=1}^{k-1} p_i^n t_{k-i}^n \quad (k \ge 2). \quad (4.1) $$

*Step 1.* Define $s_k^n = \text{Pr}(|L(\hat{G}_n)| = k)$ for $k \ge 1$. We show here that $(s_k^n)_{k \ge 1}$ is a solution to (4.1). First recall from Remark 4.2 that $\hat{G}_n$ is a.s. finite, whence $\sum_{k \ge 1} s_k^n = 1$. We introduce $G(1)$ the BCGWT planted on $\star 1$ and $\hat{G}'$ the SB-BCGWT on the right of $\star 1$ (see Figure 4.3). Clearly, for $k \ge 2$, we can write

$$ \{|L(\hat{G}_n)| = k\} = \cup_{i=1}^k \{|L(G_n(1))| = i, |L(\hat{G}'_n)| = k - i\} \cap A_n, $$

where $A_n$ is the event that there is no pruning in the edge $(\star, \star 1)$, where $G_n(1) = P_n(G(1))$ and $\hat{G}'_n = P_n(\hat{G}')$. Note that $G(1)$ and $\hat{G}'$ are independent and pruned independently, so that $G_n(1)$ and $\hat{G}'_n$ are independent. We know by Proposition 4.3 that $|L(G_n(1))|$ is $(p_k^n)_{k \ge 1}$-distributed. Furthermore, it is clear that $|L(\hat{G}'_n)|$ has the same law as $|L(\hat{G}_n)|$. Finally, exactly as in the proof of Proposition 4.3, one may check that

$$ \Pr[A_n | G_n(1), \hat{G}'_n] = \frac{2}{(|L(G_n(1))| + |L(\hat{G}'_n)| - 1)/n + 2}. $$

As a conclusion, there holds

$$ s_k^n = \sum_{i=1}^{k-1} p_i^n s_{k-i}^n \frac{2}{(i + k - i - 1)/n + 2}, $$

whence $[(k-1)/n + 2]s_k^n = 2 \sum_{i=1}^{k-1} p_i^n s_{k-i}^n$ as desired.
Step 2. Next, we show that $(kc_k^n)_{k \ge 1}$ also solves (4.1). Recall that $\sum_{k \ge 1} kc_k^n = 1$ due to (2.1). For $k \ge 2$, using (2.2),

$$[(k-1)/n + 2]kc_k^n = \sum_{i=1}^{k-1} \frac{(i + (k-i))c_i^n c_{k-i}^n}{m_0(c^n)} = 2 \sum_{i=1}^{k-1} \frac{(k-i)c_i^n c_{k-i}^n}{m_0(c^n)},$$

so that $[(k-1)/n + 2]kc_k^n = 2 \sum_{i=1}^{k-1} p_i^n (k-i) c_{k-i}^n$.

Step 3. To conclude the proof, it remains to prove that (4.1) has at most one solution. For each $x > 0$, there is obviously a unique solution $(u_k^n(x))_{k \ge 1}$ to

$$u_1^n(x) = x, \quad [(k-1)/n + 2]u_k^n(x) = 2 \sum_{i=1}^{k-1} p_i^n u_{k-i}^n(x) \quad (k \ge 2).$$

One immediately checks recursively that for each $k \ge 1$, $x \mapsto u_k^n(x)$ is increasing. Hence, there is at most one value $x_n > 0$ such that $\sum_{k \ge 1} u_k^n(x_n) = 1$. But any solution $(t_k^n)_{k \ge 1}$ to (4.1) has to satisfy $(t_k^n)_{k \ge 1} = (u_k^n(t_1^n))_{k \ge 1}$ and thus must be equal to $(u_k^n(x_n))_{k \ge 1}$. $\square$

## References

M. Aizenman and T.A. Bak. Convergence to equilibrium in a system of reacting polymers. *Comm. Math. Phys.* **65** (3), 203–230 (1979). MR530150.

D. Aldous. The continuum random tree. I. *Ann. Probab.* **19** (1), 1–28 (1991a). MR1085326.

D. Aldous. The continuum random tree. II. An overview. In *Stochastic analysis (Durham, 1990)*, volume 167 of *London Math. Soc. Lecture Note Ser.*, pages 23–70. Cambridge Univ. Press, Cambridge (1991b). MR1166406.

D. Aldous. The continuum random tree. III. *Ann. Probab.* **21** (1), 248–289 (1993). MR1207226.

G.E. Andrews, R. Askey and R. Roy. *Special functions*, volume 71 of *Encyclopedia of Mathematics and its Applications*. Cambridge University Press, Cambridge (1999). ISBN 0-521-62321-9; 0-521-78988-5. MR1688958.

P. Bak, C. Tang and K. Wiesenfeld. Self-organized criticality: an explanation of 1/f noise. *Phys. Rev. Letters* **59**, 381–384 (1987). DOI: 10.1103/PhysRevLett.59.381.

J. van den Berg and A. A. Járai. On the asymptotic density in a one-dimensional self-organized critical forest-fire model. *Comm. Math. Phys.* **253** (3), 633–644 (2005). MR2116731.

X. Bressaud and N. Fournier. On the invariant distribution of a one-dimensional avalanche process. *Ann. Probab.* **37** (1), 48–77 (2009). MR2489159.

X. Bressaud and N. Fournier. Asymptotics of one-dimensional forest fire processes. *Ann. Probab.* **38** (5), 1783–1816 (2010). MR2722786.

X. Bressaud and N. Fournier. *One-dimensional general forest fire processes*, volume 132 of *Mémoires de la Société Mathématique de France. Nouvelle Série*. Soc. Math. Fr. (2013). ISBN 978-2-85629-765-0. MR3135077.

X. Bressaud and N. Fournier. Mean-field forest-fire models and pruning of random trees. *Mém. Soc. Math. Fr.* page To appear (2014). arXiv: 1201.6645.

R. Brouwer and J. Pennanen. The cluster size distribution for a forest-fire process on $\mathbb{Z}$. *Electron. J. Probab.* **11**, no. 43, 1133–1143 (electronic) (2006). MR2268540.
J. Carr. Asymptotic behaviour of solutions to the coagulation-fragmentation equations. I. The strong fragmentation case. *Proc. Roy. Soc. Edinburgh Sect. A* **121** (3-4), 231-244 (1992). MR1179817.

D. A. Darling. On the supremum of a certain Gaussian process. *Ann. Probab.* **11** (3), 803–806 (1983). MR704564.

M. Deaconu and E. Tanré. Smoluchowski's coagulation equation: probabilistic interpretation of solutions for constant, additive and multiplicative kernels. *Ann. Scuola Norm. Sup. Pisa Cl. Sci.* **(4)** **29** (3), 549–579 (2000). MR1817709.

B. Drossel and F. Schwabl. Self-organized critical forest-fire model. *Phys. Rev. Lett.* **69**, 1629–1632 (1992). DOI: 10.1103/PhysRevLett.69.1629.

P.B. DubovSKIİ and I.W. Stewart. Trend to equilibrium for the coagulation-fragmentation equation. *Math. Methods Appl. Sci.* **19** (10), 761–772 (1996). MR1396951.

S. Finch. Bessel function zeroes (2003). Unpublished, available on http://algo.inria.fr/bsolve/.

S. Finch. Airy function zeroes (2004). Unpublished, available on http://algo.inria.fr/bsolve/.

N. Fournier and S. Mischler. Exponential trend to equilibrium for discrete coagulation equations with strong fragmentation and without a balance condition. *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.* **460** (2049), 2477–2486 (2004). MR2078689.

G. Grimmett. *Percolation*. Springer-Verlag, New York (1989). ISBN 0-387-96843-1. MR995460.

C.L. Henley. Self-organized percolation: a simpler model. *Bull. Am. Phys. Soc.* **34**, 838 (1989).

S. Janson. Brownian excursion area, Wright's constants in graph enumeration, and other Brownian areas. *Probab. Surv.* **4**, 80–145 (2007). MR2318402.

G. Louchard. The Brownian excursion area: a numerical analysis. *Comput. Math. Appl.* **10** (6), 413–417 (1985) (1984). MR783514.

R. Lyons, R. Pemantle and Y. Peres. Conceptual proofs of $L \log L$ criteria for mean behavior of branching processes. *Ann. Probab.* **23** (3), 1125–1138 (1995). MR1349164.

B. Ráth and B. Tóth. Erdős-Rényi random graphs + forest fires = self-organized criticality. *Electron. J. Probab.* **14**, no. 45, 1290–1327 (2009). MR2511285.

D. Revuz and M. Yor. *Continuous martingales and Brownian motion*, volume 293 of *Grundlehren der Mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, third edition (1999). ISBN 3-540-64325-7. MR1725357.

S. Richard and E.W. Weisstein. Catalan number (2014). From MathWorld, A Wolfram Web Resource, http://mathworld.wolfram.com/CatalanNumber.html.

P. Whittle. *Systems in stochastic equilibrium*. Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley & Sons, Ltd., Chichester (1986). ISBN 0-471-90887-8. MR850012.