The Discrete Unbounded Coagulation-Fragmentation Equation with Growth, Decay and Sedimentation

In this paper we study the discrete coagulation--fragmentation models with growth, decay and sedimentation. We demonstrate the existence and uniqueness of classical global solutions provided the linear processes are sufficiently strong. This paper extends several previous results both by considering a more general model and and also signnificantly weakening the assumptions. Theoretical conclusions are supported by numerical simulations.


Introduction
Coagulation refers to the aggregation of smaller clusters of particles to form larger ones. Some terms that are being used interchangeably with coagulation are aggregation and clustering. The first mathematical model to study such processes was proposed by Marian von Smoluchowski, who in [26,27] introduced and analysed the following system of equations (1.1) The system describes the so-called discrete coagulation, where it is assumed that any cluster consists of a finite number of monomers; that is, building blocks of minimal size (or, interchangeably, mass), taken to be equal to 1. The number of clusters of size i, called i-clusters, at any time t ≥ 0 is given by f i (t) and k i,j , i, j ≥ 1 are the coagulation rates, i.e. the rates at which the clusters of mass i and j join each other to form a cluster of mass i + j. The first term on the right-hand side, called the gain term, describes the rate of the emergence of i-clusters by coagulation of j and j − i clusters with j < i, while the second term, called the loss term, gives the rate of removal of i-clusters due to the coalescence with other ones. The factor 1 2 ensures that double counting due to symmetry is avoided.
The Smoluchowski equations describe an irreversible process. In reality coagulation is almost always coupled with a fragmentation process in which clusters split into smaller ones. The first model including fragmentation is due to Becker and Döring [7] who, however, only considered the situation in which a monomer could join or leave a cluster according to the scheme coagulation process, A comprehensive analysis of the Becker-Döring model can be found in Wattis, [28]. As far as the coagulation is concerned, the Becker-Döring model is a simplification of the Smoluchowski equation and, alleviating this shortcoming, Blatz and Tobolsky formulated a full fragmentation-coagulation model in [10]. As noted in [13], the Becker-Döring model can describe an early stage of the fragmentation-coagulation process, called the nucleation stage, when the monomers interact to build bigger clusters but still make up the majority of the ensemble. We note that there is a parallel continuous theory of fragmentation-coagulation processes in which it is assumed that the size of a particle can be any positive number. We are not concerned with such models here and the interested reader is referred to the recent monograph, [6].
Further research on fragmentation-coagulation models have led to extensions that include other internal or external processes such as diffusion or transport of clusters in space, or their decay or growth, see e.g. [11,14,29]. In particular, in applications to life sciences the clusters consist of living organisms and can change their size not only due to the coalescence or splitting, but also due to internal demographic processes such as death or birth of organisms inside, see [23,24,20,21,22]. Also, in some fields, notably in the phytoplankton dynamics, the removal of whole clusters due to their sedimentation is an important process that is responsible for rapid clearance of the organic material from the surface of the sea. The removal of clusters of suspended solid particles from a mixture is also important in water treatment, biofuel production, or beer fermentation. In all these applications the size distribution of the clusters is a crucial parameter controlling the efficacy of the process, [1,21,22]. Thus, models coupling the fragmentation, coagulation, birth, death and removal processes are relevant in many applications and hence in this paper we focus on analysing the following comprehensive system, gives the numbers f i of clusters of mass i, and, to shorten notation, we adopted the convention that g 0 = f 0 = 0. The nonnegative coefficients g i , d i and s i , i ≥ 1, control the growth, the decay and the sedimentation processes, respectively. The fragmentation rates are given by a i , while b i,j is the average number of i-mers produced after the breakup of a j-mer, with j ≥ i.
describe the rate of change of the number of particles due to, respectively, the birth and death/decay process. The form of these operators can be obtained as in the standard birth-and-death Markov process, e.g. [9], assuming that only one birth or death event can occur in a cluster of cells in a short period of time so that an i-cluster only may become an i + 1, or an i + 1-cluster. If we set g i = d i = s i = 0, i ≥ 1, then we arrive at the classical mass-conserving coagulation-fragmentation equation.
Since clusters can only fragment into smaller pieces, we have We also assume that all clusters that are not monomers undergo fragmentation; that is, a i > 0 for i ≥ 2. Since the fragmentation process only consists in the rearrangement of the total mass into clusters, it must be conservative and hence we require The main aim of this paper is to prove the existence of global classical solutions to (1.2) and provide a working numerical scheme for solving it. Thanks to recent results showing that the linear part of the problem generates an analytic semigroup, [5], in this paper we significantly extended well-posedness results existing in the literature, see e.g. [2,15], by considering more general models, removing many constraints on the coefficients of the problem and proving all results for classical solutions, and for weak solutions considered by most earlier works. The paper is organized as follows. In Section 2 we recall the main results of [5] concerning the analysis of the linear part of the problem and introduce relevant tools from the interpolation theory. Section 3 contains the proof of the global well-posedness of the problem. The idea of the analysis is classical but due to numerous technicalities specific to the problem at hand, as well as because some interim estimates are used in Section 4, we decided to provide an outline of the proofs. Finally, in Section 4 we construct finite dimensional truncations of (1.2), prove the convergence of their solutions to the solutions of (1.2) as the dimension of the truncation goes to infinity, and use the obtained results to provide rigorous numerical simulations.
Acknowledgements. The research was supported by the NRF grants N00317 and N102275, and the National Science Centre, Poland, grant 2017/25/B /ST1/00051.

The linear part
The linear part of the model (1.2) is discussed in details in [5]. In what follows we briefly mention the key results obtained there that are pertinent to the analysis of the complete nonlinear model.
In the space we consider the operators (T p , D(T p )), (G p , D(G p )), (D p , D(D p )) and (B p , D(B p )) defined by where θ i := a i + g i + d i + s i , i ≥ 1, and g 0 = a 1 = 0. Further, we denote Then the following holds (see [5] for the details): Proof. This result for p ≥ p 0 was proved in [5,Theorem 2]. Here we show that if (2.2) holds for some p 0 > 1, then it also holds for all p ∈ (1, p 0 ] and thus the argument of the proof of [5, Theorem 2] applies for all p > 1. We let It is easy to verify that for any i ≥ 2 and p > 1 0 < φ i (p) < 1. This indicates, in particular, that condition (2.2) is equivalent to the existence of constants α > 0 and β > 0 such that so that each quantity φ i (p) is monotone increasing and strictly concave on (1, ∞). The monotonicity ensures that if (2.3) is satisfied for some p 0 , then it is satisfied for and p > p 0 . On the other hand, since φ i (1) = 0, i ≥ 2, the concavity implies that for p ∈ (1, p 0 ] Hence (2.2) holds for all p > 1 and, as in [5, Theorem 2], we conclude that for each p > 1, the sum (Y p , D(Y p )) = (T p + G p + D p + B p , D(T p )) generates a positive analytic C 0 -semigroup {S p (t)} t≥0 in X p .
We mention that if the sedimentation is sufficiently strong, the generation result extends to p = 1. Indeed, in the same way as in [4,5] one can show that the assumption

Global well-posedness
In this section, we provide a well-posedness analysis of the complete semilinear model (1.2). We assume that all the conditions of Theorem 2.1 are satisfied. In addition, we impose the following bound on the coefficients of the coagulation kernel The analysis proceeds in a number of simple but technical steps. For the readers convenience the proofs of main results are broken into a sequence of short independent statements.

Local analysis
The analysis presented below is fairly standard. We convert (1.2) into an equivalent Volterra type integral equation and then employ a variant of the classical Picard-Lindelöf iterations to obtain local mild solutions. Then, with some additional work it is not difficult to verify that the mild solutions are in fact classical. The calculations are similar to that of e.g. [25, Section 6.3] or [2] but, as some intermediate estimates are needed for calculations in Section 4, we provide an outline of the proofs.
and γ = (1 + ω p + 2κ)(1 + c 0,p f 0 p,α ). As noted above, {S γ,p,α (t)} t≥0 is substochastic in X p and for all t ∈ [0, T ] and some fixed T > 0, classical solutions of (1. 2) satisfy We demonstrate that the integral equation (3.3) is locally solvable. (b) The map F γ,α : X p,α → X p is bounded and locally Lipschitz continuous provided (3.1) holds. The argument here is the same as in the proof of [2, and We use estimates (3.4) and (3.5) to show that the nonlinear map where we used the elementary inequality (3.4) and our definition of γ. Furthermore, with the aid of (3.5), (3.6), in the same manner as above we have for f, g ∈ B r (f 0 ), (d) To complete the proof we note that the maps F 1 and The last inequality indicates that B r (f 0 ) + is invariant under the action of the map M and hence the local mild solution f is non-negative.
To proceed further, we make use of the following modification of the Gronwall inequality, sometimes called the singular Gronwall inequality, see e.g. [12,Lemma 8.8.1]. Since wee need some specific aspects of it, we shall provide an elementary proof.
To simplify the notation, in the calculations below, we employ symbol c to denote a positive constant whose particular value is irrelevant. Proof. By virtue of (3.3), we have First we infer, by (2.6a) and (2.6c), For J 2 and J 3 , in the same manner as in part (c) of Lemma 3.1, we obtain Combining the estimates, we get Hence, the bound (3.16), with a constant c > 0 that depends on α, T and the initial data f 0 only, follows directly from Lemma 3.2 with γ = 1 − α. Proof. First we prove the differentiability of f in X p for t > 0. From (3.3), We observe that, by the analyticity, S γ,p,α (t)f 0 ∈ D(T p ) for t > 0, so that (3.17) The strong continuity of {S γ,p,α (t)} t≥0 , the continuity of f (see Lemma 3.1) and estimates (3.4), (3.5) combined together, show that in X p lim h→0 To find lim h→0 I 3 , we first show that where c > 0 depends on t > 0, α, T , κ, constants that appear in (2.6) and the initial data f 0 . Thus, Combining all our calculations, we conclude that df dt ∈ X p for any t > 0 and the continuity of each of the above limits shows that f ∈ C 1 ((0, T ), X p ) is a classical solution. The same calculations demonstrate also that f ∈ C((0, T ), X p,1 ).
Remark 3.5. The calculations presented above, in particular (3.17) and the last but one inequality in (3.18) Since the graph norm · p + Y p · p and · p,1 are equivalent in X p ∩ D(T p ) (as D(T p ) = D(Y p ) and both operators are closed), it follows that (3.19) for f 0 ∈ X p,α and hence f L 1 ([0,T ],Xp,1) < ∞. The last fact is crucial for the numerical analysis presented in Section 4.

Global non-negative solutions
Below we show that classical solutions of (1.2) emanating from non-negative initial data are globally defined. Our analysis requires the following elementary observation.
Lemma 3.6. Assume that f 0 ∈ X + p,α and for some ω 1 Then, under the assumptions of Theorem 3.4, the local solution satisfies Proof. Since f ∈ X + p,α , we know that every term of (1.2) is separately welldefined for t ∈ (0, T (f 0 )) (as the solution takes values in D(X p,1 )) and differentiable in X p,1 , and hence in X 1 . Thus and (3.21) follows from the standard Gronwall inequality.
Two remarks are in place here. First, in the case of pure fragmentationcoagulation models (s i = g i = d i = 0, i ≥ 1) or in the absence of growth (g i = 0, i ≥ 1), we have ω 1 ≤ 0. Second, even in the absence of sedimentation the bound (3.21) still holds provided there is a reasonable balance between the growth and the death processes. Proof. (a) To begin, we observe that for any f ∈ X + p,α we have where, by (2.2), c p and β p are positive constants that ony depend on the coefficients of (1.2) and p (in fact one can take c p to be any positive constant smaller than lim inf i→∞ ai θi with an absolute constant c 2 > 0, where we used the estimate [2, Eqn. (5.21)] for the weight. By (3.22) and the non-negativity of the local classical solution f (t), for t ∈ (0, T ) we obtain On the other hand, again by (3.1), we have while the variation of constant formula and the analiticity of the semigroup {S p (t)} t≥0 (see estimates (2.6)) imply (3.22b) We use estimates (3.22) to demonstrate that the non-negative local classical solutions cannot blow up in a finite time. For technical reasons, we separately consider two cases, 1 < p ≤ 2 and 2 < p < ∞.
Proceeding as in the proof of Lemma 3.2, we conclude that where C p,α ( f 0 p,α ) > 0 depends on the coefficients of the model (1.2), parameter 1 < p ≤ 2 and the norm f 0 p,α of the initial data, while the exponent Ω p,α > 0 is completely controlled by the parameter 1 < p ≤ 2 and the coefficients of (1.2) only. Hence, the case 1 < p ≤ 2 is settled.
(c) When 2 ≤ p < ∞, we use Hölder's inequality with the exponent q = p := p p−1 > 1. Since 0 < α < 1, we have p q (qα − 1) ≤ α, consequently and, by Young's inequality, Similar procedure yields also Hence, using (3.22a), we obtain where γ p > 0 only depends on p > 1 and the parameters of the model (1.2). From part (b) and the continuity of the embedding X 1,α ⊂ X 2,α , we have hence f (t) 1,α grows at most exponentially. Hence, the classical Gronwall inequality yields f (t) p ≤ β α,p e ω p t f 0 p also for 2 < p < ∞, where constants β α,p , ω p > 0 depend on p and the parameters of (1.2) only. As in part (b) of the proof, the last estimate, together with the inequality (3.22b), yields the exponential bound (3.23) for 2 < p < ∞. We conclude that for any p > 1, the norm f (t) p,α of the local solution f emanating from a non-negative initial datum cannot blow-up in a finite time. Hence, any such solution is defined globally. Remark 3.9. As we mentioned in Introduction, Theorem 3.7 significantly extends global solvability results obtained earlier in the context of pure coagulationfragmentation model (see [2]), where the existence of global solutions is established under much more restrictive assumptions that θ i = a i ≤ ci s , i ≥ 1, for some constants c, s > 0 and the exponent α of (3.1) satisfies 0 < αs ≤ 1.

The Truncated Problem
In numerical simulations, we approximate the original infinite dimensional system (1.2) by the following finite dimensional counterpart: (4.1) The quadratic penalty term ensures that the discrete coagulation process is conservative -this property is important when dealing with pure coagulationfragmentation models. Let P N : X p → R N and I N : R N → X p denote the projector from X p onto R N and the embedding from R N into X p , respectively. Below, we shall show that if u (N ) is the solution of the truncated problem (4.1) with the initial condition u  is non-negative, (4.2) holds for any fixed T > 0. Finally, if for some q > p−1, q ≥ 0 we have f 0 ∈ X + q+1,α and lim N →∞ I N u Proof. (a) System (4.1) is an ODE with a smooth vector field, hence it is locally solvable for any N > 0. Let where for each i ∈ N, (δ ij ) ∞ j=1 is the Kronecker delta concentrated at i. We see that the linear part of the truncated equation (4.1) acts on the elements of the finite dimensional subspace I N (R N ) ⊂ D(T p ) according to the formula Since the operator G N is non-negative and bounded, direct application of the variation of constant formula implies that the semigroup {S N (t)} t≥0 generated by (Y N , D(T p )) satisfies so that all estimates involving {S N (t)} t≥0 are uniform in N > 0. Hence, the analysis of Theorems 3.4 applies, i.e. for some T > 0 (that, in general, depends on p > 1, the initial condition and the coefficients of the problem) inclusion (4.2) holds and the respective norms are bounded independently of N .
Assuming that the initial datum u (N ) 0 is non-negative, we proceed as in Theorem 3.7 to show that the inclusion (4.2) holds for any fixed T > 0 uniformly in N . Hence, the first two claims of Theorem 4.1 are settled.
(b) To prove the last claim, we derive the equation governing evolution of the numerical error e (N ) (t) : where, for a given f and u (N ) , H N (t)e (N ) is linear in e (N ) and In what follows we will use two inequalities based on the properties of the function [0, a] x → φ(x) := x r (a − x) r , a > 2, r > 0. Clearly, φ is symmetric, nonnegative with φ(0) = φ(a) = 0 and has a single maximum at x = a/2. Thus, for x ∈ [1, a − 1] we have φ(x) ≥ (a − 1) r . In particular, for q ≥ 0 and a = N + 1 we have where the inequality for q = 0 is trivial, and for p ≥ 1, using (4.4) and j ≤ N Then, by (3.1), (3.4), the fact that f is globally defined by Theorem 3.7, and (4.5), H N e (N ) satisfies where (4.5) was used to get Similarly, using (4.4) to estimate E 1 N f , we have where all generic constantsc > 0 are uniform in N > 0. The last four bounds, combined with the variation of constants formula, with a constantc > 0 that does not depend on the truncation parameter N > 0. Further, by (4.6) and (3.19), we have where, as before,c > 0 is independent of N > 0. Thus, using (

Simulations
Below, we provide several numerical illustrations to the theory developed above. In our simulations, we make use of the following two fragmentation kernels: Example 2. In our second example, we employ the fragmentation kernel (4.7b) with σ = 10 −1 and the coagulation kernel (4.8b) with k 2 = 5 · 10 −3 and k 3 = 1. Note that k i,j = O(i 2 + j 2 ) and, in view of (3.1), we let δ = 2.5. The remaining set of parameters is identical to those used in Example 1.
In the settings described above, the growth rate of the quantities k i,j is superlinear. Hence, the pure coagulation models lead to a formation of a massive particle outside the system (the so called gelation phenomenon, see [28] and references therein). In addition, the moment conditions, proposed in [2] in context of the discrete pure coagulation-fragmentation models, are also not satisfied. Nevertheless, the example fells in the scope of Theorem 3.7 and, as predicted by the theory, the numerical solution demonstrates qualitative features simi- 2) with the coagulation kernel (4.8b) and the fragmentation kernel (4.7b): number of clusters u n (t) (top left); distribution of cluster masses nu n (t) (top right); the total number of particles (middle left); the total mass (middle right) and the higher order moments (bottom). lar to those observed in Example 1, see Fig. 2. The total mass is preserved (i.e. no shattering and/or gelation occur) and after a short transition stage the numerical trajectory settles near a stationary particles/mass distribution. chosen to be the same as in Examples 1 and 2.

The growth-decay-sedimentation-fragmentation-coagulation scenario
As demonstrated by Fig. 3, in the presence of the transport processes the qualitative dynamics of the model (1.2) changes (compare Fig. 3 with Fig. 1  and 2). The death and the sedimentation processes dominate and yield a slow decay in each of the moments u p , p = 0, 1, 2, 3 as time increases.
Example 4. To provide a further illustration of the effect of transport processes on the dynamics of (1.2), we repeat the computations but with the fragmentation and the coagulation kernels from Example 2. To ensure global solvability of the model, we let g = d = s = a = 1, β = γ = 0 and α = δ = 2.5.
With this settings, the birth and the fragmentation terms dominate and we expect the total mass of the ensemble to grow. As shown in Fig. 4, this is indeed the case for t close to zero. However, as time goes on, the contributions of the growth and the decay/sedimentation processes compensate each other and the numerical solution settles near an equilibrium state. The example demonstrates certain degree of flexibility of model (1.2). A proper interplay between the fragmentation and the transport components of the equation allows for simulation of a wide range of realistic scenarios arising within coupled transport-fragmentation-coagulation systems.  Figure 5: Evolution of the decay-sedimentation-coagulation-fragmentation model (1.2) with the coagulation kernel (4.8a) and the fragmentation kernel (4.7a): number of clusters u n (t) (top left); distribution of cluster masses nu n (t) (top right); the total number of particles (bottom left) and the total mass (bottom right).

The no-growth scenario
Our last two examples demonstrate behaviour of (1.2) in the absence of growth, i.e. when g = 0 and with sufficently strong sedimentation. In this settings, the model is globally well posed in X 1 , provided (2.4) and (3.1) are satisfied.
Example 5. We let g = 0, d = s = a = 1, γ = δ = 1 and β = 0. The fragmentation and the coagulation kernels and all other parameters are the same as in Example 1.
The results of simulations are shown in Fig. 5. The strong sedimentation (see condition (2.4)) describing the death of clusters, prevents uncontrolled mass absorption by the clusters of large sizes. The top right diagrams in Fig. 5 demonstrate that the bulk mass of the system remains concentrated in clusters of moderate size. As time goes on, both processes lead to a steady decay in the total mass of the system. Example 6. In our last example, we make use of the fragmentation and the coagulation kernels from Examples 2 and 4. Further, we set g = 0, d = s = a = 1, γ = δ = 2.5 and β = 0.
As mention earlier, the growth rate of the quantities k i,j is superlinear and one expects gelation in context of pure coagulation models. Nevertheless, in  Figure 6: Evolution of the decay-sedimentation-coagulation-fragmentation model (1.2) with the coagulation kernel (4.8b) and the fragmentation kernel (4.7b): number of clusters u n (t) (top left); distribution of cluster masses nu n (t) (top right); the total number of particles (bottom left) and the total mass (bottom right).
complete agreement with the theory, the simulations show (see the evolution of the clusters masses in the top right diagrams of Fig. 6) that in the presence of a sufficiently strong decay-sedimentation process the latter scenario is impossible, and the solution remains bounded in X 1 settings (see the bottom right diagram in Fig 6). It is worth to mention that in this example the mechanism preventing gelation is connected with the strong sedimentation, in contrast to Examples 2 and 4 where the central role is played by the strong fragmentation.

Conclusion
In the paper, we considered the discrete coagulation-fragmentation models with growth, decay and sedimentation. The analysis presented in Section 3 shows that, irrespective of the coagulation rates, the model is always globally well posed, provided the fragmentation (in the case of p > 1), or the sedimentation (for p = 1) dominate. This is in contrast to pure coagulation models, see e.g. [28]) but confirms earlier results obtained in a more restricted setting in the discrete, [2,15], and continuous, [18], cases. Theoretical conclusions are completely supported by the numerical simulations presented in Section 4.