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

HAL
open science

# Asymptotic behaviour of a one-dimensional avalanche model through a particular stochastic process

Madalina Deaconu, Oana Lupaşcu-Stamate

► To cite this version:

Madalina Deaconu, Oana Lupaşcu-Stamate. Asymptotic behaviour of a one-dimensional avalanche model through a particular stochastic process. 2022. hal-03947249

HAL Id: hal-03947249

https://inria.hal.science/hal-03947249

Preprint submitted on 19 Jan 2023

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.
# Asymptotic behaviour of a one-dimensional avalanche model through a particular stochastic process

Madalina Deaconu,* and Oana Lupascu-Stamate†‡

## Abstract

In this paper we develop the study of a binary coagulation-fragmentation equation which describes the avalanches phenomena. We construct first an adapted stochastic process and obtain its behaviour to the equilibrium. Our model is based on self-organized critical (SOC) systems and in particular on a simple sand pile model introduced in Bressaud and Fournier, [BrFo 09]. Furthermore, we define a stochastic differential equation for this process and propose a numerical method in order to approximate the solution. The key point of our work is a new interpretation of the avalanches phenomena by handling stochastic differential equations with jumps and the analysis of the invariant behaviour of the stochastic process.

**Keywords:** Fragmentation and coagulation equation, stochastic differential equation, self-organized critical system, avalanches modelling, numerical methods.

2010 MSC: 34A34, 60H15, 65N75, 74S60, 74G15.

## 1 Introduction

In this paper we study a binary coagulation-fragmentation equation which can be used for modelling of avalanches phenomena. Our results open new directions in this topic. In particular the enhancement of stochastic models and the construction of numerical algorithms are opening new perspectives for associated real complex problems.

Self-organized critical (SOC) systems are very rich and important tools that describe complex phenomena. The concept was introduced first by Bak, Tang and Wiensefeld [BaTaWi 87] in the framework of the sand pile model. SOC can be interpreted as the property of dynamical systems to organize their microscopic behaviour. They are in general simple models that describe spatial and temporal randomness in phenomena showing long range correlations. They can be used to study various real phenomena as: avalanches, earthquakes, fires, financial market crashes, rock landslips, etc.

For short, the concept introduced in Bak, Tang and Wiensefeld [BaTaWi 87] for the sand pile is the one of a dynamical system expanding toward a critical state and conducing to a rupture and

<footnote>*Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France.</footnote>

<footnote>†Gheorghe Mihoc-Caius Iacob Institute of Mathematical Statistics and Applied Mathematics of the Romanian Academy, Calea 13 Septembrie 13, RO-050711 Bucharest, Romania</footnote>

<footnote>‡Research Group of the Project PN-III-P4-PCE-2021-0921 Simion Stoilow Institute of Mathematics of the Romanian Academy, Bucharest, Romania</footnote>

<footnote>Email: Madalina.Deaconu@inria.fr, oana.lupascu_stamate@yahoo.com</footnote>
giving catastrophic events. In this evolution, when the behaviour becomes more and more probable
while approaching the critical state, then the dynamical system reaches an equilibrium which is close
to this critical state.

The SOC state can be also observed in the classical model of a given pile where the grains can
be added randomly or at some fixed time. As a consequence of these events, the system can became
unstable, and thus can be interpreted and modelled as an avalanche or a disaster.

By using these properties, in the paper [OlFeCh 92] the authors constructed a description for
the earthquakes. Based also on the SOC behaviour there are also many works for the forest fire
evolution, see [He 89], [DrSc 92], [BrFo 10].

Starting with a simple sand model, in [BrFo 09] an avalanche model was introduced in the
following sense: if a grain falls on a position that is not empty, then the entire connected component
including this position becomes free.

We will consider in our approach this class of SOC models. In [BrFo 09], the authors investigated
the equilibrium properties for the avalanche process, from theoretical and numerical point of view.
First, it is proved that the avalanche process has an unique invariant distribution with an adapted
rate of convergence.

Secondly, a deterministic coagulation-fragmentation equation related to the avalanche process,
is considered. For the pure coagulation, the origin of this equation goes back to the important work
of Smoluchowski in 1916, see [Sm 16].

Let us describe the evolution of the system under consideration: for $i \in \mathbb{N}^*$ it is written down an
infinite system of ordinary differential equations satisfied by $c(i, t)$, the concentration of the clusters
of mass $i$ at time $t$ (a cluster of mass $i$ is formed by $i$ identical particles, each particle having a mass
1). Specific coagulation and fragmentation kernels drive the dynamic of the process.

Numerical results show that the equilibrium state of this coagulation-fragmentation system is a
good approximation of the invariant distribution of the avalanche process [BrFo 09].

Partial results on the trend to equilibrium for binary coagulation-fragmentation models without
balance condition are presented in [FoMi 04].

For the discrete Smoluchowski's coagulation equation with some solvable coagulation kernels,
in [MePe 03] it is proved the uniform convergence of densities to the self-similar solution with
exponential tails under adapted moment hypotheses.

In the paper [EsMiRi 05], under the self-similarity property for the coagulation-fragmentation
equation, is obtained the existence of a stationary solution for any given mass and it is shown
that the solution of the fragmentation equation behaves, as $t \to \infty$, as the unique self-similar
solution of the same mass. In [LaPe09] for constant coefficients in the linear growth-fragmentation
equation, it is proven that the dynamic converges to the steady state with an exponential rate.
The authors in [BaJoSh 19] showed that the long-term behaviour of any solution of the discrete
growth-decay-fragmentation equation, is given by a scalar exponential function multiplied by a
vector independent on the initial conditions.

Recent results [BeDeLu 15, BeDeLu 16, BeDeLu-St 19] are considering the link between the frag-
mentation equation and branching processes. The approach of the fragmentation with probabilistic
tools (SDE, branching processes) is also developed. In particular we are getting inspired in our
numerical approach by the algorithm constructed for the simulation of these stochastic processes in
[BeDeLu-St 19].

In this paper, we introduce a binary coagulation-fragmentation equation for an avalanche phe-
nomenon based on self-organized critical (SOC) system. A SOC system is defined by an interacting
particle system taking values in $\{0, 1\}^\mathbb{Z}$, see [BrFo 09] for more details. We introduce specific coag-
ulation and fragmentation kernels related to the particular behaviour observed in the avalanche phenomenon. Moreover, using an implicit formulation of the steady state of the coagulation-fragmentation system, we study its equilibrium properties (see Proposition 3.1 from Section 3).

Further, we develop a probabilistic approach for our particular coagulation-fragmentation equation, by constructing an adapted stochastic differential equation with pure jumps (see equation (4.4) (SDECF)). Based on this probabilistic interpretation, we introduce a numerical approximation for the avalanche process, Algorithm A1 from Section 5, and we establish a convergence result in Theorem 5.1. Then, comparing this with the values for the steady state of the initial deterministic coagulation-fragmentation equation, we construct an approximation (see Table 1) in Section 5.

In this framework, our results are stronger and finer than the existing results in the literature on this topic. A particular analysis is needed for the fragmentation part and also for the construction of a good balance between coagulation and fragmentation evolution.

The paper is organized as follows. After the introduction part, in Section 2 we construct an avalanche model, following the framework in [BrFo 09], described by an interacting particle system taking values in $\{0, 1\}^Z$.

For this model, in Section 3 we construct an alternative binary coagulation-fragmentation equation in discrete time, by introducing an appropriate binary phenomenon with specific coagulation and fragmentation kernels, and then we describe the link with the model introduced in Section 2. Our model mixes both the classical coagulation-fragmentation equation and the particular behaviour observed in the avalanche model. The particularity of our fragmentation procedure is that a particle $i$ could be only split into two particles one of mass 1 and the other one of mass $i-1$. This is similar, in some sense, to the avalanche phenomenon from Section 2, see the description of the connection between the coagulation-fragmentation equation considered and the avalanche model. At the end of the Section 3 we deduce the equilibrium properties of the avalanche process by using an implicit formulation of the steady state of the coagulation-fragmentation system.

In particular, in Section 4 we construct a new and more general numerical approximation for the avalanche model introduced in Section 2, by considering a probabilistic approach for the specific coagulation-fragmentation equation from Section 3. Our model is based on a stochastic differential equation with pure jumps. We consider that the conservation of mass holds under some hypothesis that ensure existence and uniqueness of the solution (cf. [FouGi 03] and [Jou 03]).

The interpretation of the evolution of the stochastic differential equation is that the mass of a typical particle in the system is obtained either by adding, at some random Poissonian times, driven by the coagulation kernel $K$, the mass of another typical particle (if coagulation occurs) ; either by splitting, at some random Poissonian times, driven by the fragmentation kernel $F$, this mass into two smaller masses. This allows us to construct a new numerical approximation for the solutions of such stochastic differential equations.

In Section 5 we construct a recursive algorithm to approximate the solution of the stochastic differential equation of coagulation-fragmentation introduced in Section 4, based on its probabilistic interpretation and we establish a convergence result. We also give some numerical results for this algorithm for the avalanche model. Then, comparing them with the values for the steady state of the initial deterministic coagulation-fragmentation equation, we obtain a good approximation. Finally, we give a convergence result of the algorithm based on its recursive description.
## 2 An avalanche model

We construct an avalanche model which is given by an interacting particle system with values in $\{0, 1\}^Z$. This construction is based on the model introduced in [BrFo 09].

Let $\Gamma = ((\Gamma_t(i))_{t \ge 0})_{i \in \mathbb{Z}}$ be an independent family of Poisson processes with rate 1. In this model, we assume that on each site $i \in \mathbb{Z}$, the snowflakes are falling following the Poisson process $(\Gamma_t(i))_{t \ge 0}$ and the birth flakes follows Poisson processes with rate 1.

For a site $i \in \mathbb{Z}$ we say that it is empty if $i = 0$ and occupied when $i = 1$. If a flake falls on a vacant site $i$ of $\mathbb{Z}$, i.e. $i = 0$, this site becomes occupied, $i = 1$. If a flake falls on an occupied site $i \in \mathbb{Z}$, an avalanche starts, that means that the whole connected component of occupied sites around $i$ becomes vacant.

Let us mention that we can also construct more general models in which the killing flakes follow Poisson processes with rate greater than 1, while the birth flakes follow Poisson processes with rate 1. In this case, forest fire models can be considered. These models were studied for example in [BrFo 10, BrFo 13, BrFo 14].

We can now define the avalanche process and denote it by $((y_t(i))_{t \ge 0})_{i \in \mathbb{Z}}$. The definition for $t \ge 0$ and $i \in \mathbb{Z}$, is as follows:

* $y_t(i) = 1$ if the site $i$ is occupied at time $t$, and,

* $y_t(i) = 0$ if the site $i$ is vacant at time $t$.

The state space of this process is the following:

$$ (2.1) \qquad E := \{\mathbf{y} \in \{0, 1\}^Z \text{ such that } \liminf_{i \to -\infty} \mathbf{y}(i) = \liminf_{i \to \infty} \mathbf{y}(i) = 0\}. $$

Notice that the condition from the definition of $E$ insures that the states of the process have no infinite connected component with occupied sites for an initial condition $y_0$ belonging to $E$. However, as explained in [BrFo 09], if $y_0 \in \{0, 1\}^Z \setminus E$, then there exists $t_0 > 0$, such that $y_t \in E$ for any $t \ge t_0$.

For any initial condition in $E$ and for almost any family of Poisson processes $\Gamma$ as before, the avalanche process $\mathbf{y} = (y_t)_{t \ge 0}$ exists, it is unique, and it has the strong Markov property, see for example [BrFo 09] and [Lig 05].

We shall consider this framework for our study.

## 3 A related binary coagulation-fragmentation model and numerical algorithm

Let us now introduce a binary coagulation-fragmentation equation for the avalanche model considered in Section 2. We describe in this section the link between the two models.

The discrete coagulation-fragmentation equation is an infinite system of non-linear ordinary differential equations describing the time evolution of a system of particles in which each particle is characterised by its mass. The equation gives the evolution in time of $c(t, i)$ which is the concentration
tion of particles of mass $i$ at time $t$ in the system. The equation takes the form

$$ (3.1) \qquad \left\{ \begin{aligned} \frac{\partial}{\partial t} c(i,t) &= \frac{1}{2} \sum_{j=1}^{i-1} (K(i-j,j)c(i-j,t)c(j,t) - F(i-j,j)c(i,t)) \\ & \quad - \sum_{j \in \mathbb{N}^*} (K(i,j)c(j,t)c(i,t) - F(i,j)c(i+j,t)) \\ c(i,0) &= c_0(i), \end{aligned} \right. $$

for $(i,t) \in \mathbb{N}^* \times \mathbb{R}_+$. In this model the particles can either coagulate to form bigger particles or fragment into smaller ones. We consider here only binary phenomena.

The coagulation and fragmentation kernels $K$ and $F$ are nonnegative symmetric functions defined on $\mathbb{N}^* \times \mathbb{N}^*$ with positive real values. Here $K(i,j)$ represents the coagulation rate of two particles with masses $i$ and $j$ to form a cluster with mass $i+j$, while $F(i,j)$ is the fragmentation rate of a cluster $i+j$ that breaks into two particles with masses $i$ and $j$. We refer to [Sm 16, DaCo 93, Je 97, Jou 03] for more details on the discrete coagulation-fragmentation process.

The first fragmentation term in (3.1) counts the particles of size $i$ which disappears after splitting into two smaller particles $i-j$ and $j$, while the second fragmentation term is counting the creation of particles of size $i$, after the breakage of particles of larger size $i+j$, into two particles $i$ and $j$. The first coagulation term represents the particles of size $i$ resulting after the coagulation of two particles $i-j$ and $j$. The second coagulation term counts the particles $i$ that disappear after coagulating with an other particle $j$.

From the description of the model we can expect solutions to have the property of the conservation of mass, that means:

$$ \sum_{i \ge 1} ic(i,t) = \sum_{i \ge 1} ic_0(i) \quad \text{for all } t \ge 0, $$

where we suppose that the initial distribution has finite mass: $\sum_{i \ge 1} ic_0(i) < \infty$. However, explicit solutions in the pure coagulation case ($F \equiv 0$) for the kernel $K(i,j) = ij$ show that this is not true. In this case the mass is not conserved over the time. In particular at a time called gelification time the mass decreases, and this is interpreted as the formation of a cluster with infinite mass, see [Al 99] and the references cited therein.

In [Je 97] the problem of gelation is treated by using probabilistic tools by approximating the solutions of the coagulation-fragmentation equations with a sequence of finite state Markov chains in $l_2$ space. In the pure fragmentation case ($K \equiv 0$) it is possible to construct solutions with increasing mass, see [BaCa 90].

In our case, we consider solutions of coagulation-fragmentation equation with non-increasing total mass, in the sense of [Jou 03, DaCo 93], as follows.

For $T \in (0, \infty]$, we say that the map $t \in [0, T) \to c(\cdot, t) \in \{c : \mathbb{N}^* \to \mathbb{R}_+, \sum_{i \ge 1} ic(i) \le \sum_{i \ge 1} ic_0(i)\}$ solves the coagulation-fragmentation equation (3.1) on $[0, T)$ if for all $i \in \mathbb{N}^*$ and for all $t \in [0, T)$, the applications $u \to \sum_{j \in \mathbb{N}^*} K(i,j)c(j,u)$ and $u \to \sum_{j \in \mathbb{N}^*} F(i,j)c(i+j,u)$ are integrable on $(0,t)$ and the following equation holds:

$$ (3.2) \qquad c(i,t) = c_0(i) + \int_0^t \frac{1}{2} \sum_{j=1}^{i-1} (K(i-j,j)c(i-j,u)c(j,u) - F(i-j,j)c(i,u)) \\ \qquad \quad - \sum_{j \in \mathbb{N}^*} (K(i,j)c(j,u)c(i,u) - F(i,j)c(i+j,u))du. $$
Now, we construct an alternative coagulation-fragmentation equation by introducing an appropriate binary phenomenon with specific coagulation and fragmentation kernels. This coagulation-fragmentation equation is connected to (3.1) and to the avalanche model introduced in Section 2. Let us consider the following coagulation and fragmentation kernels: for all $i, j \in \mathbb{N}^*$ we take $K(i, j) = 2$ and

$$ (3.3) \qquad F(i, j) = \begin{cases} i + j - 1, & \text{if } i = 1 \text{ or } j = 1, \\ 0, & \text{if } i \neq 1 \text{ and } j \neq 1. \end{cases} $$

The particularity of the fragmentation kernel is that a particle $i$ could be only split into two particles one of mass 1 and the other one of mass $i-1$. This is similar to the avalanche phenomenon studied by [BrFo 09].

We consider in the front of the coagulation term an acceleration factor $1/m_0(t)$ depending on the total number of particles per unit of length, i.e. $m_0(t) = \sum_{i \ge 1} c(i, t)$, see [BrFo 09] for more details. The equation (3.1) becomes

$$ (3.4) \qquad \left\{ \begin{aligned} \frac{\partial}{\partial t} c(i,t) &= \frac{1}{m_0(t)} \sum_{j=1}^{i-1} c(i-j,t) c(j,t) - (i+1)c(i,t) + ic(i+1,t) \\ c(i,0) &= c_0(i). \end{aligned} \right. $$

We explain in this part the link between the equation (3.4) and the avalanche model described in Section 2. We refer also to [BrFo 09] for alternative presentation of this. Moreover, we introduce a mean-field model related to the problem inspired by the presentation in [BrFo 09]. The basis is the same but we propose here a particular mathematical model tailor-made for the avalanche that we study.

Let $(\mathbf{y}_t)_{t \ge 0}$ be the avalanche process with state space $E$, described in Section 2. In particular this choice will insure that the construction of the coagulation-fragmentation process associated to (3.4) will preserve the total mass.

For a state of the avalanche process $\mathbf{y}_t \in E$, we say that two neighbour edges $(i-1, i)$ and $(i, i+1)$ belong to the same particle (cluster) if $\mathbf{y}_t(i) = 1$ (the site $i$ is occupied at time $t$). So, a particle of mass $i$ contains (covers) $i$ edges and $i-1$ occupied sites. In particular, $(i, j)$ belongs to a particle with mass 1 if and only if $\mathbf{y}_t(i) = \mathbf{y}_t(j) = 0$. We assume that each edge has a mass equal to 1.

For a configuration $\mathbf{y} \in E$, and for $i \in \mathbb{N}$, we assume that there exists the average number of particles with mass $i$ per unit of length denoted by $c(i, \mathbf{y})$.

We suppose in a simplified framework that we have independence between the successive masses of particles in $\mathbf{y}_t$ at each time $t \ge 0$. We admit thus that

$$ (3.5) \qquad c(i,t) := c(i, \mathbf{y}_t) = \frac{1}{i} \mathbb{P}[\text{the edge } (0,1) \text{ belongs to a particle with mass } i \text{ in } \mathbf{y}_t]. $$

Then $(c(\cdot, t))_{t \ge 0} = (c(i, t))_{i \ge 0, t \ge 0}$ verifies the corresponding mass conservation property $\sum_{i \ge 0} ic(i, t) = 1$ for all $t \ge 0$. Moreover, $(c(\cdot, t))_{t \ge 0} = (c(i, t))_{i \ge 0, t \ge 0}$ is solution of the following infinite system of
differential equations, which is similar to the coagulation-fragmentation equation (3.4):

$$ (3.6) \quad \left\{ \begin{aligned} \frac{\partial}{\partial t} c(1,t) &= -2c(1,t) + \sum_{i \ge 1} ic(i+1,t), \\ \frac{\partial}{\partial t} c(i,t) &= -2c(i,t) - (i-1)c(i,t) + \frac{1}{m_0(t)} \sum_{j=1}^{i-1} c(i-j,t)c(j,t) + ic(i+1,t) \quad \text{for all } i \ge 2. \end{aligned} \right. $$

Our model mixes both the classical coagulation-fragmentation equation and the particular behaviour observed in the avalanche model (for the particles of mass 1). Thus the case $i = 1$ has to be treated separately.

The first equation in (3.6) counts the particles with mass 1 (the isolated edge). The first term on the right says that such a particle disappears when a flake falls on an isolated edge and it merges (coagulates) at a constant rate 1 with its two neighbours. This process corresponds to the constant coagulation rate $K(i, j) = 2$, for all $i, j$. The second term is counting the phenomenon of fragmentation that means a particle with mass 1 appears if a flake falls on a particle of mass $i+1$ ($i+1$ edges and $i$ sites) and an avalanche starts: it breaks into two particles of masses 1 and $i$. That event corresponds to the fragmentation rate $F(i, 1)$ defined by (3.3).

In the second equation of (3.6), the first term on the right side represents the particles of mass $i \ge 2$ which disappear after becoming larger with rate 2, when a flake falls on one of its extremities. The second term counts the particles of mass $i$ which disappear when a flake falls on a particle with mass $i$ and an avalanche occurs with the rate $F(i - j, j) = i - 1$ defined by (3.3), that means the particle of mass $i$ breaks into two smaller particles of masses $i - j$ and $j$. The third term counts the particles of mass $i$ which appear when a flake falls on one extremity of a particle of mass $j$, which is the extremity of a particle of mass $i - j$ and they merge (coagulate with the constant rate 2). Finally, the last term counts the particles of mass $i$ which result after the splitting of a larger particle $i+1$ into two smaller particles (of masses $i$ and 1), according to the fragmentation kernel $F(i, 1)$.

The objective is now to develop the equilibrium properties of the avalanche process $y_t$ by using the steady state of the system (3.6) and a stochastic model for the equation (3.4).

Let us start by proving an implicit formulation of the steady state of the system (3.6).

**Proposition 3.1.** The system of equations (3.6) admits a unique steady state denoted by $c = (c(i))_{i \ge 1}$, that is:

$$ (i) \quad c(i) \ge 0 \text{ for each } i \ge 1 $$

$$ (ii) \quad \sum_{i \ge 1} ic(i) = 1 $$

$$ (iii) \quad c = (c(i))_{i \ge 1} \text{ satisfies} $$

$$ (3.7) \quad 2c(1) = \sum_{i \ge 1} ic(i+1), $$

$$ (3.8) \quad (i+1)c(i) = ic(i+1) + \frac{1}{M_c} \sum_{j=1}^{i-1} c(j)c(i-j) \quad \text{for all } i \ge 2, $$

where $M_c := \sum_{i \ge 1} c(i)$ is supposed uniquely defined and finite.
The steady state is given implicitly by (3.8). We can approximate numerically $M_c \simeq 0,459134$ and then the steady state of the system.

Proof. Since $M_c = \sum_{i \ge 1} c(i)$, using condition (ii) and the equation (3.7), we get that

$$ (3.9) \qquad c(1) = \frac{1 - M_c}{2}. $$

Clearly, $M_c \ge 0$. From condition (i), we deduce $M_c \in [0, 1]$.

For $i \ge 2$, $c(i)$ is uniquely determined from the recurrence (3.8) and condition (3.9) as a rational function on $M_c$ whose denominator is a power of $M_c$. More precisely, $c(2i - 1)$ and $c(2i)$ are two polynomials, where the numerator has the degree $i$, while the degree of the denominator is $i - 1$ due to the term $\frac{1}{M_c}$ from (3.8).

By computation $c(2) = 1 - 2M_c$. Since $c(2) \ge 0$, adding to the previous restriction we get $M_c \in [0, \frac{1}{2}]$.

Continuing with $c(3) = \frac{-25M_c^2 + 14M_c - 1}{8M_c} \ge 0$, we obtain $M_c \in [0.084, 0.47]$.

Continuing with $c(4) = \frac{-29M_c^2 + 20M_c - 3}{6M_c} \ge 0$, we get $M_c \in [0.220, 0.469]$.

We used the Singular software ([DGPS 22]) to compute the real roots of the equation $c(i) = 0$ up to $i = 60$. The condition (i) yields $M_c \in [0.459134112, 0.459134206]$.

We conclude that $M_c \simeq 0,459134$. Using the equation (3.8) we can numerically compute $c = (c(i))_{i \ge 1}$. For more numerical results see the Section 5. By hypothesis $M_c$ is uniquely determined, thus the steady state of the equation (3.6) is also uniquely determined. $\square$

# 4 Probabilistic approach for the avalanche model

We construct a new and more general numerical approximation for the avalanche model introduced in Section 2 by considering a probabilistic approach for the specific Smoluchowski coagulation-fragmentation equation (3.4) from Section 3.

We can suppose that, without losing the generality, that $\sum_{i \ge 1} ic_0(i) = 1$ (plays the role of a normalisation constant). Thus $p(i) := ic_0(i)$, for every $i \in \mathbb{N}^*$, is a probability distribution on $\mathbb{N}^*$.

Further, we introduce a probabilistic approach for the discrete coagulation-fragmentation equation (3.4) that can be used to model the avalanche. We use some ideas here from the model introduced and developed in [FouGi 03] for continuous masses.

For more generality, we develop in this part the continuous formulation of the coagulation-fragmentation equation which can be adapted it in the discrete case. The probabilistic interpretation of the general coagulation-fragmentation equation is naturally constructed by using the conservation of mass property. Let us denote by $\mu(t, x) = xc(t, x)$ for every $t \ge 0$ and $x \in \mathbb{N}^*$. Then the coagulation-fragmentation equation with continuous masses writes:

$$ (4.1) \qquad \begin{cases} \frac{\partial}{\partial t}c(t, x) = \frac{1}{2} \int_0^x c(t, y)c(t, x - y)K(y, x - y)dx - c(t, y) \int_0^{+\infty} c(t, y)K(x, y)dy \\ -\frac{1}{2} \int_0^x c(t, x)F(y, x - y)dy + \int_x^{+\infty} F(x, y - x)c(t, y)dy \\ c(0, x) = c_0(x). \end{cases} $$
For the discrete masses situation, with a particular choice of the distribution, we can prove that the solution of the continuous coagulation-fragmentation equation (4.1) is also solution of the discrete coagulation-equation (3.1) (for more details we explain below this property).

We introduce the following definition:

**Definition 4.1.** Let $T \le +\infty$ be fixed. We say that the family of nonnegative measures $(c(t, \cdot))_{t \in [0,T)}$ on $\mathbb{R}_+$ is a weak solution to the coagulation-fragmentation equation if

* for all $t \in [0, T)$, $\int_0^{+\infty} xc(t, dx) = 1$

* for all $t \in [0, T)$, $\sup_{s \in [0,t]} \int_0^{+\infty} \int_0^{+\infty} xyK(x, y)c(s, dx)c(s, dy) < \infty$

* for all $t \in [0, T)$ and all test function $\varphi : \mathbb{R}_+ \to \mathbb{R}$ such that $\varphi(x)/x \in C_b^1(\mathbb{R}_+)$

$$ \int_0^{+\infty} \varphi(x)c(t, dx) = \int_0^{+\infty} \varphi(x)c(0, dx) + \frac{1}{2} \int_0^t \int_0^{+\infty} \int_0^{+\infty} [\varphi(x+y) - \varphi(x) - \varphi(y)]K(x, y)c(s, dx)c(s, dy)ds + \frac{1}{2} \int_0^t \int_0^{+\infty} \int_0^x [\varphi(x-y) + \varphi(y) - \varphi(x)]F(y, x-y)c(s, dx)dyds. $$

We consider here that the conservation of mass holds and, as stated above, this assumption allows to write in a compact form the discrete and the continuous mass situations and give a probabilistic interpretation of the coagulation-fragmentation equation. The conservation of mass writes in the discrete mass situation

$$ (4.2) \qquad \sum_{i=1}^{+\infty} ic(i, t) = \sum_{i=1}^{+\infty} ic_0(i) = 1 $$

while in the continuous case it writes

$$ (4.3) \qquad \int_0^{+\infty} xc(x, t)dt = \int_0^{+\infty} xc_0(x)dx = 1. $$

We observe that $Q_t(dx) = \sum_{i=1}^{+\infty} ic(i, t)\delta_i(dx)$ is a probability measure on $\mathbb{N}$ and also $Q_t(dx) = xc(t, x)dx$ is a probability measure on $\mathbb{R}_+$ for each $t$. Thus the discrete and continuous mass situations lead to the same equation in terms of $Q_t$.

We aim now to give a stochastic interpretation of the solution as a stochastic differential equation with pure jumps and need to introduce some hypothesis on the coefficients of the equation.

The form of our equation can be used on a $G = \mathbb{N}^*$ or $G = \mathbb{R}_+$. We keep the notation $G$ in order to have a general framework.

Our study is concentrated on the discrete case thus for us $G = \mathbb{N}^*$. We recall a standard hypothesis that ensure the existence and the uniqueness of the solution (deterministic case). We will use this solution to construct the solution of the SDE.
We state first the following hypothesis :

**Hypothesis 1**

1. The initial condition $c_0(dx)$ is a nonnegative finite measure such that $\int_G xc_0(dx) = 1$.

2. The coagulation kernel $K: G \times G \mapsto \mathbb{R}_+$ is a continuous symmetric map. There exists a constant $C$ such that for all $x, y \in G$

$$K(x, y) \leq C(1 + x + y).$$

There exists a continuous nonnegative function $\phi: G \to [1, \infty)$ such that $x \to \phi(x)/x$ is nonincreasing on $G$, and for all $x, y \in G$,

$$0 \leq K(y, x) = K(x, y) \leq \phi(x)\phi(y).$$

We also have, $\int_G (x^2 + \phi^2(x))c_0(dx) < \infty$ and either $x \to \phi^2(x)/x$ is nonincreasing on $G$ or for all $x, y \in G$, $K(x, y) \leq \phi(x) + \phi(y)$.

3. The fragmentation kernel $F: G \times G \mapsto \mathbb{R}_+$ is a continuous symmetric map. The function $\Psi: G \mapsto \mathbb{R}_+$ defined by

$$\Psi(0) = 0, \Psi(x) = \frac{1}{x} \int_{0}^{x} y(x - y)F(y, x - y)dy, \text{ for } x > 0,$$

is continuous. There exist constants $p \in \mathbb{N}^*$ and $C > 0$ such that $\Psi(x) \leq C(1 + x^p)$. Furthermore, for all $\varepsilon \in (0, 1)$

$$\lim_{n \to \infty} \sup_{x \in [\varepsilon, \frac{1}{\varepsilon}]} \frac{1}{x} \int_{0}^{x} y(x - y)F(y, x - y)1_{F(y, x - y) \geq n}dy = 0.$$

We introduce the notion of solution for the stochastic differential equation of coagulation-fragmentation (SDECF) in the following definition.

**Definition 4.2.** A stochastic process $(X_t)_{t \geq 0}$ is a solution of the (SDECF) if there exists a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \geq 0}, \mathbb{P})$ such that

1. $X_0$ is a $xc_0(dx)$ distributed random variable,

2. $(X_t)_{t \geq 0}$ is a càdlàg G-valued $(\mathcal{F}_t)_{t \geq 0}$ adapted process,

3. There exist two independent Poisson measures $N(ds, dy, dz)$ and $M(ds, dy, dz)$ adapted to $(\mathcal{F}_t)_{t \geq 0}$ on $[0, +\infty) \times G \times [0, +\infty)$ with intensity measures $dsQ_s(dy)dz$, and $dsdydz$ respectively, where for $s \geq 0$, $Q_s$ is the distribution of $X_s$, such that a.s. for all $t \geq 0$

$$(4.4) \qquad \begin{aligned} X_t = X_0 &+ \int_0^t \int_0^{+\infty} \int_0^{+\infty} y1_{\{z \leq \frac{K(X_{s-}, y)}{y}\}} 1_{\{y < \infty\}} N(ds, dy, dz) \\ &- \int_0^t \int_0^{+\infty} \int_0^{+\infty} y1_{\{y \in (0, X_{s-})\}} 1_{\{z \leq \frac{X_{s-} - y}{X_{s-}} F(y, X_{s-} - y)\}} M(ds, dy, dz). \end{aligned}$$
Under the Hypothesis 1 there exists a weak solution to the (SDECF), see the result in [FouGi 03].
In the discrete case similar results are developed in [Jou 03].

This formulation in terms of the stochastic differential equation is of great interest as we can con-
struct a numerical algorithm by using the associated jump Markov process. We give here the general
framework for the existence of the solution of the stochastic differential equation of coagulation-
fragmentation (4.4).

However in our situation we are in the discrete case and the hypothesis are not satisfied by our
particular fragmentation kernel. To overcome this difficulty, as we have simple particular forms for
the kernels, we will construct directly the solution, by developing a recurrence procedure (see section
5).

**Intuitive idea** : The main idea of the evolution of the stochastic differential equation is to consider a particle system with some given initial masses with the well-chosen distribution. The mass of a typical particle in the system is obtained either by adding, at some random Poissonian times, driven by the coagulation kernel $K$, the mass of another typical particle (if coagulation occurs) ; either by splitting, at some random Poissonian times, driven by the fragmentation kernel, this mass into two smaller masses. This remark allows us to construct a recursive algorithm to approximate the stochastic process $X_t$.

Based on this interpretation, we develop the following Algorithm to simulate the solution of the
stochastic differential equation of coagulation-fragmentation (4.4).

5 Algorithm

We develop in this section the numerical method constructed for the approximation of the solution
of the coagulation-fragmentation equation (4.4) and prove a convergence result of the algorithm. A
result concerning only the coagulation phenomenon was considered in Proposition 3.2 in [FouGi 04].
Our result is stronger and more difficult to obtain as it needs particular attention to the fragmentation
part and the construction of a good balance between coagulation and fragmentation evolution.

We denote in the following, for $a > 0$, by $\mathcal{E}(a)$ the exponential distribution with parameter $a$.

Algorithm A1 - Simulation of the coagulation-fragmentation SDE

Step 0: Sample $X_0^1, \dots, X_0^n$ i.i.d.r.v., according to $Q_0(dx) = xc_0(dx)$.

Fix the final time $T$ and set $T_0 = 0$.

Step p.1: Sample a random variable $y^i \sim \mathcal{U}([0, X_{T_{p-1}}^i])$ for all $i \in \{1, \dots, n\}$.

Compute $m_{p,f}^i = \frac{X_{T_{p-1}}^i - y^i}{X_{T_{p-1}}^i} F(y^i, X_{T_{p-1}}^i - y^i)$ for all $i \in \{1, \dots, n\}$.

Sample a random variable $S_{p,f}^i \sim \mathcal{E}(m_{p,f}^i)$ for all $i \in \{1, \dots, n\}$.

$S_{p,f} = \min_i \{S_{p,f}^i\}, \quad i_f = \arg \min_i \{S_{p,f}^i\}.$

Step p.2: Compute $m_{p,c} = \sup_{i,j} \frac{K(X_{T_{p-1}}^i, X_{T_{p-1}}^j)}{X_{T_{p-1}}^j}$ for all $i,j \in \{1, \dots, n\}$.

Set $i_c = i, j_c = j$ the couple for which the sup is realised.

Sample a random variable $S_{p,c} \sim \mathcal{E}(nm_{p-1,c})$.
**Step p.3:** Set $S_p = \min\{S_{p,f}, S_{p,c}\}$.

**Step p.3f:** If $S_p = S_{p,f}$ then $X_{T_{p-1}}^{i_f}$ could fragment.

Sample a random variable $u \sim \mathcal{U}([0, 1])$.

If $u \le m_{p,f}^{i_f}$, then fragmentation occurs and set

$$X_{T_p}^{i_f} = X_{T_{p-1}}^{i_f} - y^{i_f}.$$

$$X_{T_p}^k = X_{T_{p-1}}^k \quad \text{for all } k \neq i_f.$$

Else set $X_{T_p}^k = X_{T_{p-1}}^k$, for all $k \in \{1, \dots, n\}$.

**Step p.3c:** If $S_p = S_{p,c}$ then a coagulation can occur at time $S_p$.

Sample $u \sim \mathcal{U}([0, m_{p,c}])$.

If $u \le \frac{K(X_{T_{p-1}}^{i_c}, X_{T_{p-1}}^{j_c})}{X_{T_{p-1}}^{j_c}}$, then the coagulation occurs and set

$$X_{T_p}^{i_c} = X_{T_{p-1}}^{i_c} + X_{T_{p-1}}^{j_c}$$

$$X_{T_p}^k = X_{T_{p-1}}^k \quad \text{for all } k \neq i_c.$$

Else set $X_{T_p}^k = X_{T_{p-1}}^k$ for all $k \in \{1, \dots, n\}$.

**Step p.4:** Set $T_p = T_{p-1} + S_p$ and for all $k \in \{1, \dots, n\}$ and $t \in [T_{p-1}, T_p)$,

$$X_t^k = X_{T_{p-1}}^k.$$

**Stop:** When $T_p > T$, set, for all $t \in [T_{p-1}, T]$, $X_t^k = X_{T_{p-1}}^k$ for all $k \in \{1, \dots, n\}$.

**Outcome:** The approximated particles mass at time $T$, $X_{T_{p-1}}^k$ for all $k \in \{1, \dots, n\}$.

## Numerical results.

In Table 1 we give the values of $c_i, i = 1, \dots, 5$, the concentration of the particles of mass $i$ in the system, obtained with the Algorithm A1 for the fragmentation kernel (3.3) and the coagulation kernel $K(i, j) = \frac{2}{M_c}$ with $M_c = 0.459134$ (from Proposition 3.1), $N = 15$, $T = 100$ and Monte Carlo parameter $M = 10^4$. Then, we compare them with the values obtained in Proposition 3.1.

<table><caption>Table 1: Numerical approximations of concentration of particles $c_i, i = 1, \dots, 5$</caption><thead><tr><th>c<sub>i</sub></th><th>Algorithm A1</th><th>Value given by Proposition 3.1</th></tr></thead><tbody><tr><td>c<sub>1</sub></td><td>0.2707</td><td>0.270433</td></tr><tr><td>c<sub>2</sub></td><td>0.1003</td><td>0.081732</td></tr><tr><td>c<sub>3</sub></td><td>0.0209</td><td>0.042954</td></tr><tr><td>c<sub>4</sub></td><td>0.0050</td><td>0.025178</td></tr><tr><td>c<sub>5</sub></td><td>0.0013</td><td>0.015186</td></tr></tbody></table>
In this particular situation we conclude that our algorithm gives a good approximation for the steady state of the deterministic coagulation-fragmentation equation (3.6) related to avalanches.

Moreover, using (3.8) from Proposition 3.1 and the Matlab software, we can compute more approximations for $c_6 \approx 0.009232$, $c_7 \approx 0.005625$, $c_8 \approx 0.003430$, $c_9 \approx 0.002092$, $c_{10} \approx 0.001275$, $c_{11} \approx 0.000778$, $c_{12} \approx 0.000474$, $c_{13} \approx 0.000289$, $c_{14} \approx 0.000176$, $c_{15} \approx 0.000107$, $c_{16} \approx 6.568997e-05$, $c_{17} \approx 4.006588e-05$, $c_{18} \approx 2.443714e-05$, $c_{19} \approx 1.490479e-05$, $c_{20} = 9.090790e-06$, ..., $c_{50} = 2.960125e-12$.

So, we can compute $\sum_{i=1}^{50} c_i = 0.459036$, which is a good approximation of $M_c$.

## 5.1 The convergence of the algorithm

We start this section by rewriting our Algorithm A2 in a recursive way. Consider for this the following algorithm:

**Algorithm A2 - Recursive Algorithm**

Sample $X_0^1, \dots, X_0^n$ i.i.d.r.v., according to $Q_0(dx) = xc_0(dx)$.

Fix the final time $T$ and set the initial $t=0$;

Set $x = X_0 = (X_0^1, \dots, X_0^n)$;

While $(t < T)$, do

{

Compute the fragmentation or coagulation time $S_p = \min\{S_{p,f}, S_{p,c}\}$ and the $i_f, i_c, j_c$ as in Step p.1 and Step p.2 from the Algorithm A1.

Set $X_s = x$ for all $s \in [t, (t + S_p) \wedge T]$

Set $t = t + S_p$

If $t \le T$

{

Choose $w$ uniformly in $[0, 1]$

Set $xm=\text{mass}(t,w,x)$

Set $x(i_f) =xm(i_f) - y(i_f)$, for the fragmentation case

or set $x(i_c) =xm(i_c) +xm(j_c)$, for the coagulation case.

}

}

We follow here some ideas from Proposition 3.2 in [FouGi 04], in the coagulation pure case and we define a random recursive function $\text{mass}(t,v,x)$ which gives the procedure to obtain the mass of a particle in the fragmentation / coagulation system as follows: we consider the current particle in the algorithm. At the next step in the fragmentation case we subtract from this particle another small particle (chosen randomly) with an adapted fragmentation rate, or with an adapted coagulation rate we add to this particle another particle from the system (chosen randomly), in the coagulation case.

This function is defined by the following procedure:

function $\text{mass}(t, v, z)$

{

Simulate $X_0$ a random variable distributed as $Q_0(dx)$.

Set $x = X_0$ and $s = 0$
While $s < t$ and test of coagulation ($v \le K(x, z)/z$)
or fragmentation ($v \le m_{p,f}^{if}$) (see Step p.1, Algorithm 1) valid do
$
\begin{aligned}
& \text{Compute the fragmentation or coagulation time } S_p = \min\{S_{p,f}, S_{p,c}\} \\
& \text{Set } s = s + S_p \\
& \text{If } s \le t \\
& \quad \begin{cases}
\text{Choose } w \text{ uniformly in } [0, 1] \\
\text{Set } y = \text{mass}(s, w, x) \\
\quad \text{If } S_p = S_{p,f}, \text{ a fragmentation can occur:} \\
\quad \quad x = x - y \\
\quad \text{else } S_p = S_{p,c}, \text{ a coagulation can occur:} \\
\quad \quad x = x + y
\end{cases} \\
& \quad \text{If } v \le K(x, z)/z, \text{ cf. Step p.3c} \\
& \quad \quad \text{set } \text{mass}(r, v, z) = x \text{ else} \\
& \quad \quad \text{set } \text{mass}(t, v, x) = 0 \\
& \quad \text{If } v \le m_{p,f}^{if}, \text{ cf. Step p.3c} \\
& \quad \quad \text{set } \text{mass}(t, v, z) = x \text{ else} \\
& \quad \quad \text{set } \text{mass}(t, v, z) = 0 \\
& \quad \quad p = p + 1. \\
& \quad \end{cases}
\end{aligned}
$
}

We prove next a convergence result of our algorithm, based on its recursive description. In our case, the constant kernel $K(x, y) = 2$ fullfills the conditions 1 and 2 from the Hypothesis 1.

For a measure $\nu$ and a function $f$, we denote by

$$< \nu(dx), f(x) > = \int f(x) d\nu(x).$$

We give now the convergence result of the algorithm to the solution of (4.4).

**Theorem 5.1.** Let $T < \infty$. We assume that $K$ satisfies the conditions 1 and 2 from the Hypothesis 1 and suppose that for all $x, x', y \in \mathbb{N}^*$ if $x \le x'$ then $K(x, y) \le K(x', y)$. We suppose also that $K$ and $F$ are such that there exist a solution to the (SDECF). Denote by $C_T$ the total number of times that Recursive Algorithm takes before ending.

Then

$$ (5.1) \qquad \mathbb{E}[C_T] < \infty $$

and the Recursive Algorithm ends a.s. We denote by $(X_t)_{t \in [0,T]}$ be the process constructed by the algorithm. Then $(X_t)_{t \in [0,T]}$ satisfies the equation (4.4).

*Proof.* Let $T < \infty$ be fixed.

For any $r \ge 0$, any $v \in (0, 1)$, and any $z \in \mathbb{N}^*$, we denote by $P_{r,v,z}(dx, dc)$ the law of a couple of random variables $(X_{r,v,z}, C_{r,v,z})$, with $C_{r,v,z}$ the (possibly infinite) number of times that the execution
of the recursive function mass(r,v,z) calls the function mass, and $X_{r,v,z}$ is the result of the function mass(r,v,z). By convention we set $X_{r,v,z} = 0$ if $C_{r,v,z} = \infty$. Notice that the function mass(r,v,z) gives the possibility to have a fragmentation or a coagulation procedure.

For each $r \ge 0$, we denote by $C_r$ the (possibly infinite) total number of times that Recursive Algorithm calls the function mass to obtain $X_r$. Then $C_r$ is a nondecreasing $\mathbb{N} \cup \{\infty\}$-valued process. We set $X_r = \infty$ if $C_r = \infty$.

Then, since $X_{r,v,z}$ is simulated essentially in the same way as $X_r$, in law, we have

$$ (X_{r,v,z}, C_{r,v,z}) \stackrel{(d)}{=} (X_r(1_{\{v \le [K(z,X_r)/X_r]/\lambda(z)\}} 1_{\{X_r < \infty\}} 1_{\text{coag.}} - 1_{\{X_r \in (0,z)\}} 1_{\{u \le \frac{z-X_r}{z} F(X_r,z-X_r)\}} 1_{\text{frag.}}, C_r), $$

where $\lambda(z) = m_{p,c}$ from Step p.2 in the Algorithm A1.

For each $r \in [0, T]$ we denote by $Q_r$ the law of the $\mathbb{N}^* \cup \{\infty\}$-valued random variable $X_r$.

The process $(X_r)_{r \in [0,T]}$ is now well-defined as a càdlàg, $\mathbb{N}^* \cup \{\infty\}$-valued process, and $X_0$ has the distribution $xQ_0(dx)$.

Let $D$ be the set $\mathbb{N} \times (\mathbb{N} \cup \{\infty\})$. For each $r \in [0, T]$ we have that

$$ (5.3) \qquad X_r = X_0 + \int_0^r \int_0^1 \int_D x M^c(ds, dv, d(x, c)) - \int_0^r \int_0^1 \int_D x M^f(ds, du, d(x, c)), $$

where $M^c(ds, dv, d(x, c))$ and $M^f(ds, du, d(x, c))$ are random integer-valued measures (see [JacSh 87]) on $[0, T] \times (0, 1) \times D$, with compensator $\lambda(X_{s-})dsdvP_{s,v,X_{s-}}(dx, dc)$, respectively $dsduP_{s,u,X_{s-}}(dx, dc)$. We also have

$$ (5.4) \qquad C_r = \int_0^r \int_0^1 \int_D (1+c) M^c(ds, dv, d(x, c)) + \int_0^r \int_0^1 \int_D (1+c) M^f(ds, dv, d(x, c)). $$

Further, we show that $(X_r)_{r \in [0,T]}$ satisfies (SDECF).

Using (5.2), we can rewrite (5.3) using two random integer-valued measures $O^c(ds, dv, dx)$ and $O^f(ds, du, dx)$ on $[0, T] \times (0, 1) \times (\mathbb{N}^* \cup \{\infty\})$ with compensator $\lambda(X_{s-})dsdvQ_s(dx)$, respectively $dsduQ_s(dx)$, as

$$ (5.5) \qquad X_r = X_0 + \int_0^r \int_0^1 \int_{\mathbb{N}^* \cup \{\infty\}} x 1_{\{v \le \frac{K(X_{s-},x)}{\lambda(X_{s-})x}\}} 1_{\{x < \infty\}} O^c(ds, dv, dx) \\ - \int_0^r \int_0^1 \int_{\mathbb{N}^* \cup \{\infty\}} x 1_{\{x \in (0, X_{s-})\}} 1_{\{u \le F(x, X_{s-}-x) \frac{X_{s-}-x}{X_{s-}}\}} O^f(ds, du, dx). $$

Using two Poisson measures $N(ds, dy, dv)$ and $M(ds, dy, du)$ on $[0, T] \times (\mathbb{N}^* \cup \{\infty\} \times [0, \infty))$ with compensator $\lambda(X_{s-})dsQ_s(dydv)$, respectively $dsQ_s(dydu)$, we can rewrite the above stochastic differential equation in the form of the equation (SDECF), given by (4.4).

We prove now that:

$$ (5.6) \qquad \sup_{t \in [0,T]} \mathbb{E}[X_t + \lambda(X_t)] < \infty. $$
Let us remark that for the case restricted to the coagulation equation this result was given in Lemma 2.4 from [FouGi 04]. In our case we use in particular the inequality (3.6) in the above mentioned paper.

Indeed let us denote by $X_t$ the solution to the equation (SDECF) and denote by $X_t^c$ the solution of the equation of pure coagulation with the same coagulation kernel, driven by the same Poisson measure.

Suppose that we start with the same initial condition $X_0 = X_0^c$.

For $X_t^c$ denote by $T_1^c$ the first time when a coagulation occur.

For $X_t$ we have to compare two possible times : the coagulation time $T_1^c$ (the same as before !) and the fragmentation time $T_1^f$.

* If $T_1^c < T_1^f$ then $X_t = X_t^c$ for all $t \in [0, T_1^c]$ and the analysis is going on from this new starting point.

* If $T_1^f < T_1^c$ then $X_t = X_t^f$ for all $t \in [0, T_1^f)$ and $X_{T_1^f} \le X_{T_1^c}$.

The main observation is that if we start with an initial condition $X_0 \le X_0^c$ then we can prove, step by step, by using the form of the stochastic differential equation and the hypothesis on the $K$, that at any time

$$ (5.7) \qquad \forall t \in [0, T], X_t \le X_t^c. $$

By using the Lemma 2.4 from [FouGi 04] (for the coagulation case) and (5.7) we deduce that the result is also valid for the process $X_t$. Thus, with (5.6) we deduce that, almost surely

$$ (5.8) \qquad C_T < \infty. $$

This implies that the algorithm ends almost surely in a finite number of steps.

Let us now prove that $\mathbb{E}[C_T] < \infty$.

For $r \in [0, T]$ let $n_r$ be the distribution of the random variable $C_r$ which is $\mathbb{N} \cup \{+\infty\}$ valued. By using formula (5.2) we can express (5.4) on the following form

$$ (5.9) \qquad C_r = \int_0^r \int_{\mathbb{N}} (1+c)\nu(ds, dc), $$

where $\nu$ is a random integer-valued measure of the form $\nu = \nu^c + \nu^f$, and $\nu^c$ and $\nu^f$ denote also random integer-valued measures with compensators $\lambda(X_{r-})drn_r/dc$ and $drn_r/dc$ respectively.

Let $A \in (0, +\infty)$, we have:

$$ (5.10) \qquad C_r \wedge A \le \int_0^r ds((1+c) \wedge A)\nu(ds, dc). $$

By taking expectation, we get

$$ (5.11) \qquad \begin{aligned} \mathbb{E}(C_r \wedge A) &\le \int_0^r ds \mathbb{E}(\lambda(X_s) + 1) \int_{\mathbb{N}} [(1+c) \wedge A] n_s/dc \\ &\le \int_0^r ds \mathbb{E}(\lambda(X_s) + X_s) \int_{\mathbb{N}} [(1+c) \wedge A] n_s/dc. \end{aligned} $$
Thus, by denoting $m_T = \sup_{t \in [0,T]} \mathbb{E}(\lambda(X_t) + X_t)$ we obtain by using (5.11)

$$ (5.12) \qquad \mathbb{E}(C_r \wedge A) \le m_T + m_T \int_0^r \mathbb{E}(C_T \wedge A)ds. $$

And finally, by applying Gronwall Lemma we deduce that $\mathbb{E}(C_T \wedge A) \le \gamma_T$, where $\gamma_T$ is a constant not depending on $A$. Taking the limit as $A$ goes to infinity leads to the conclusion.

This ends the proof of the convergence result.

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

Let us conclude that we have developed a new probabilistic procedure that gives an approximation of the coagulation-fragmentation equation using stochastic differential equations with jumps, for particular kernels, that can be used to approach the avalanche phenomenon and opens new perspectives.

**Acknowledgements.** For the second author this work was supported by a grant of the Ministry of Research, Innovation and Digitization, CNCS - UEFISCDI, project number PN-III-P4-PCE-2021-0921, within PNCDI III.
References

[Al 99] D. J. Aldous, Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the mean-field theory for probabilists, *Bernoulli* **5**(1) (1999), 3–48.

[BaTaWi 87] P. Bak, C. Tang, K. Wiesenfeld, Self-organized criticality: An explanation of 1/f noise, *Phys. Rev. Lett.* **59** (1987), 381–284.

[BaCa 90] J. M. Ball, J. Carr, The discrete coagulation-fragmentation equations: existence, uniqueness and density conservation, *Journal of Statistical Physics* **61** (1990), 203–234.

[BaJoSh 19] J. Banasiak, L.O.Joel, S. Shindin, Discrete growth–decay–fragmentation equation: well-posedness and long-term dynamics, *Journal of Evolution Equations* **19** (2019), 771–802.

[BeDeLu 15] L. Beznea, M. Deaconu, O. Lupaşcu, Branching processes for the fragmentation equations, *Stochastic Processes and their Applications* **125** (2015), 1861–1885.

[BeDeLu 16] L. Beznea, M. Deaconu, O. Lupaşcu, Stochastic equation of fragmentation and branching processes related to avalanches, *Journal of Statistical Physics*, **162** (2016), 824–841.

[BeDeLu-St 19] L. Beznea, M. Deaconu, O. Lupaşcu-Stamate, Numerical approach for stochastic differential equations of fragmentation, application to avalanches, *Mathematics and Computers in simulation* **160** (2019), 111–125.

[BrFo 09] X. Bressaud, N. Fournier, On the invariant distribution of a one-dimensional avalanche process, *Annals of Probability* **37** (2009), 48–77.

[BrFo 10] X. Bressaud, N. Fournier, Asymptotics of one-dimensional forest fire processes, *Annals of Probability* **38** (2010), 1783–1816.

[BrFo 14] X. Bressaud, N. Fournier, A Mean-field forest-fire model, *ALEA, Lat. Am. J. Probab. Math. Stat.* **11**(2) (2014), 589–614.

[BrFo 13] X. Bressaud, N. Fournier, Processus de feux de forêt généraux en dimension 1, *Mémoires de la SMF* (2013), DOI:10.24033/msmf.442.

[DaCo 93] F.P. Da Costa, Existence and uniqueness of density conserving solutions to the coagulation-fragmentation equations with strong fragmentation, *Journal of Mathematical Analysis and Applications* **192** (1995), 892–914.

[DGPS 22] Decker, W.; Greuel, G.-M.; Pfister, G.; Schönemann, H.: SINGULAR 4-3-0 — A computer algebra system for polynomial computations. https://www.singular.uni-kl.de (2022).

[DrSc 92] B. Drossel, F. Schwabl, Self-organized critical forest-fire model, *Phys. Ref. Lett.* **69** (1992), 1629–1632.
[EsMiRi 05] M. Escobedo, S. Mischler, M. Rodriguez Ricard, On self-similarity and stationary problem for fragmentation and coagulation models, *Ann. I. H. Poincaré* **22** (2005), 99–125.

[FouGi 03] N. Fournier, J.-S. Giet, On small particles in coagulation-fragmentation equations, *Journal of Statistical Physics* **111** (2003), 1299–1329.

[FouGi 04] N. Fournier, J.-S. Giet, Exact simulation of nonlinear coagulation processes, *Monte Carlo Methods and Applications* **10**(2) (2004), 95–106.

[FoMi 04] N. Fournier, 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** (2004), 2477–2486.

[He 89] C. L. Henley, Self-organized precolation: A simpler model, *Bull. Amer. Phys. Soc.* **34** 838 (1989).

[JacSh 87] J. Jacod and A.N Shiryaev, *Limit Theorems for Stochastic Processes*, Springer Verlag, (1987).

[Je 97] I. Jeon, Existence of gelling solutions for the coagulation-fragmentation equations, *Communications in Mathematical Physics* **194** (1998) 541–567.

[Jou 03] B. Jourdain, Nonlinear processes associated with the discrete Smoluchowski coagulation-fragmentation equation, *Markov Processes and Related Fields* **9** (2003), 103–130.

[LaPe09] P. Laurencot, B. Perthame, Exponential decay for the growth-fragmentation/cell-division equation, *Communications in Mathematical Sciences*, International Press **7**(2) (2009), 503–510.

[Lig 05] T.M. Liggett, *Interacting Particle Systems*, Springer, 2005.

[MePe 03] G. Menon, R. L. Pego, Dynamical Scaling in Smoluchowski's Coagulation Equations: Uniform Convergence, *SIAM Review* **36**(5) (2003), DOI: 10.2307/20453875.

[OlFeCh 92] Z. Olami, H. J. S. Feder, K. Christensen, Self-organized critically in a continuous, non-conservative cellular automation modeling earthquakes, *Phys. Rev. Lett.* **68** (1992), 1244–1247.

[Sm 16] M. Smoluchowski, Drei Vortrage über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen, *Physik. Zeitschr.* **17** (1916), 557–599.