Steady oscillations in aggregation-fragmentation processes

We report surprising steady oscillations in aggregation-fragmentation processes. Oscillating solutions are observed for the class of aggregation kernels $K_{i,j} = i^{\nu}j^{\mu} + j^{\nu}i^{\mu}$ homogeneous in masses $i$ and $j$ of merging clusters and fragmentation kernels, $F_{ij}=\lambda K_{ij}$, with parameter $\lambda$ quantifying the intensity of the disruptive impacts. We assume a complete decomposition (shattering) of colliding partners into monomers. We show that an assumption of a steady-state distribution of cluster sizes, compatible with governing equations, yields a power law with an exponential cutoff. This prediction agrees with simulations results when $ \theta \equiv \nu-\mu<1$. For $ \theta=\nu-\mu>1$, however, the densities exhibit an oscillatory behavior. While these oscillations decay for not very small $\lambda$, they become steady if $\theta$ is close to two and $\lambda$ is very small. Simulation results lead to a conjecture that for $ \theta<1$ the system has a stable fixed point, corresponding to the steady-state density distribution, while for any $\theta>1 $ there exists a critical value $\lambda_c(\theta )$, such that for $\lambda<\lambda_c(\theta)$, the system has an attracting limit cycle. This is rather striking for a closed system of Smoluchowski-like equations, lacking any sinks and sources of mass.


I. INTRODUCTION
Numerous phenomena in nature involve dual processes of aggregation and fragmentation [1,2]. These processes take place on vastly different length and time scales. A reversible polymerization in solutions and coagulation of colloidal particles are the classical examples of such processes occurring on the molecular scales; another peculiar example is aggregation of prions causing the Alzheimerlike diseases [3]. On the larger scales-in atmospheric processes, small airborne particles coalesce into smog droplets [4]. Aggregation is also common in systems of living organisms, from colonies of viruses [5] to schools of fish [6]. Aggregation and fragmentation processes occur in networks of different nature, including economic networks [7] and internet communities [1,8]; here forums of users nucleate, merge and split. In turbulent cascades in a fluid flow [9], vortices may merge forming larger ones or decomposing into smaller vortices. The distribution of particles size in planetary rings is also determined by a steady balance achieved between two opposite processes, viz. aggregation and breakage of the particles in the rings [10][11][12][13][14].

A. Aggregation
The aggregation takes place when two clusters, comprised respectively of i and j monomers, merge upon collision thereby creating a cluster of i + j monomers (see Fig. 1); symbolically this process may be written as where K ij is the merging rate. Let n k be the concentration of clusters of size k, i.e., clusters composed of k monomers. The rate of change of n k is determined by Smoluchowski equations [1,2] The first term in the right-hand side accounts for the formation rate of k-mers from clusters of size i and j, the second term describes the loss of k-mers due to aggregation of these clusters with all other clusters; the factor 1/2 in the first term prevents from double counting of the same process (i + j → k and j + i → k).

B. Aggregation with fragmentation
Generally aggregates can suffer both spontaneous and collision fragmentation [1-3, 10, 11, 13]. In the former case a cluster breaks into smaller pieces without interactions with other aggregates [1][2][3], in the latter one the fragmentation is caused by an energetic impact between two clusters [10,11,13]. Different collision fragmentation models have been studied [10,11,13]; here we will consider a simple one of complete shattering of two colliding 2 partners into monomers. Symbolically this process (see Fig. 1

) may be written as [i] + [j]
Fij − − → [1] + [1] + . . . [1] i+j where F ij quantifies the shattering rate. Models with shattering exhibit interesting behaviors including dynamical phase transitions [15]. It has been shown [10] that more general fragmentation models with large number of fragments yield qualitatively similar size distribution provided the small-size debris strongly dominates over the large size ones [10]. As in Ref. [10], we assume that the fragmentation and aggregation kernels are proportional, as has been justified for processes in planetary rings [10]. The parameter λ in Eq. (2) characterizes the relative frequency of aggregative and shattering impacts.
Adding the fragmentation kinetics with kernel (2) into kinetic equations (1) we arrive at a separate equation for the concentration of monomers and a set of generic equations for k ≥ 2. The second term on the right-hand side of Eq. (3) accounts for the gain of monomers occurring in shattering collisions between clusters, the third term describes the gain of monomers in the shattering impacts between monomers and clusters. Equation (4) differs from (1) by an extra loss term (proportional to λ) accounting for shattering. Equations (3)-(4) describe spatially-homogeneous systems. The kernels K i,j may be obtained from the microscopic analysis of the aggregation and fragmentation processes (see e.g. [10,13,16]). In applications, K i,j are usually homogeneous functions of the masses i and j of merging clusters. Here we will investigate the kernels which are rather popular [1,2] and have been used for similar aggregating-shattering systems in Ref. [17] where a source of monomers and sink of large aggregates was present. A stationary distribution satisfying Eqs. (3)-(4) with the kernel (5) has been also addressed in [18]. In the following, we shall often use the sum and the difference of the exponents µ and ν (Without loss of generality, we choose ν ≥ µ.) The exponent β is the well-known homogeneity exponent [1,2]. The exponent θ plays an important role in the following; it has been called a non-locality exponent in [17,18]. We always limit ourselves to the non-gelling case β < 1. The restrictions ν ≤ 1 and µ ≤ 1 are needed to avoid instantaneous gelation (see e.g. [19][20][21][22][23]). The exponent θ can exceed 1, and never-ending oscillations are actually observed in a 'non-local' regime θ > 1.
In the special case of ν = −µ = a, the kernel reads with 0 ≤ a ≤ 1. This kernel is known as a generalized Brownian kernel [24]. In what follows we will analyze both (5) and (7), often starting with the latter which is more tractable. The restriction a ≤ 1 is needed to avoid instantaneous gelation (aggregation equations with a > 1 are ill-defined [19][20][21][22][23]25]). The solutions of the aggregation-fragmentation equations should also satisfy the natural physical requirement n k (t) ≥ 0, and mass conservation: Analytical time-dependent solutions to Eqs. (3)-(4), have been obtained only for the simplest case of a constant kernel [10]. The steady-state solutions have been found for several other models, such as irreversible aggregation model with a monomer source [26], an aggregationfragmentation model with kernels K i,j = (ij) µ and F i,j = λK i,j [10], and for an open aggregation-fragmentation system with a source of monomers and sink of large clusters [17] for the kernels of the form (5) and for closed systems in [18]. An open aggregating system with the same coagulation kernel (5), driven by input of monomers along with the removal of large clusters has been studied in [27]. Steady oscillations were numerically found in this system with a finite number of aggregate species [27]. For a closed system comprised of monomers, dimers, trimers and exited monomers, stable oscillations have been also reported [28]. Similarly, steady chemical oscillations may occur in a dimerization model (see, e.g., [29]). In the present study we consider closed systems undergoing aggregation and fragmentation processes, with the kinetic rates given by Eqs. (5) and (2), that lack any source or sink of monomers and clusters. Naively, one expects that such closed systems with two opposite processes will relax to a steady-state where a balance between aggregation and shattering is established. This scenario is indeed realized for θ < 1, or a < 1/2 in the case of the kernels (7). Unexpectedly, for θ → 2 (or a → 1) and small values of λ we observe never-ending oscillations of the concentrations. This effect has been found numerically for the one-parameter family of kernels (7) and reported in our recent study [30]. Here we present a more detailed analysis of the aggregating and shattering systems, both numerical and theoretical, and we investigate a more general two-parameter family of kernels (5). We also provide a qualitative theory of the stable oscillations which sheds some light on the mechanism of this surprising phenomenon.
In what follows, we will concentrate on systems with time-independent coefficients K ij and λ, and conserved total mass (total number of elementary units). This is a generic model that describes systems of very different nature. The elements comprising a system range from grains or molecules to living organisms or economic agents. The interaction forces, that determine the kinetic rates, may be of very different nature as well. These may range from true molecular or mechanical forces to fictitious "social forces" [31] based on informational exchange. Therefore, strictly speaking, the term "closed" here literally means a lack of sinks and sources of system elements. At the same time, the exchange of energy, chemicals, nutrients, and information is implied.This is needed to sustain elements of a system and keep the rate coefficients steady. On the level of social agents or a living organism this implies an interaction with the surrounding social or natural environment. On the level of molecular or macroscopic particles the interaction with a thermostat, or the presence of some other source of energy, is assumed. For instance, in aggregation-fragmentation processes in polymer or colloidal solutions, there exists energy exchange with the solvent. This maintains constant temperature, although energy is released in aggregation processes and consumed in fragmentation processes. Similarly the surrounding molecular gas plays a role of thermostat in atmospheric processes [32] and for dust clouds [33,34]. Another important mechanism of energy supply is viscous heating, which arises in planetary rings [14]. In this case the orbital motion of rings' particles yields a sheared flow of viscous granular fluid, which generates heat [14]. The energy supply keeps the kinetic energy of aggregates steady, and the rate coefficients constant.
Systems with true molecular or mechanical forces between elements is an important subclass of systems with aggregation and fragmentation. As it follows from the discussion above, such systems, with constant rate coefficients, are not thermo-dynamically closed. (Note that the notion "thermodynamics" is meaningful only for these systems.) Hence an interesting question ariseswhether persistent concentration oscillation sexist in thermodynamically closed systems? We perform a microscopic analysis, which resulted in a positive answer: Never-ending oscillations do emerge in thermodynamically closed systems, although the oscillation period permanently increases.
The rest of the paper is organized as follows. In the next Sec. II, we present simulation results obtained with the use of fast solvers of Smoluchowski-type equations. In Sec. III we discuss steady-state distributions using the methods outlined in Ref. [30] and applied to Brownian kernels µ = −ν = a. We also present a qualitative theory explaining the mechanism leading to never-ending oscillations. In Sec. IV we summarize our findings.

II. NUMERICAL RESULTS
Kinetic equations (3)-(4) form a set of infinitely many nonlinear coupled ordinary differential equations (ODE), which is a severe numerical challenge. For standard Smoluchowski equations, that is when fragmentation is absent, the average size of aggregates grows indefinitely imposing a time limit to model these processes. Fragmentation precludes the formation of very large clusters (in most cases and certainly in our case when fragmentation and aggregation kernels are proportional). This allows to model the aggregating-and-shattering systems with a finite number of equations N eq , which is dictated by the requested accuracy. In Ref. [30] we present estimates that relate the number of equations and the simulation accuracy; in practice we use such number of equations, that a further increase of N eq does not impact the results for the concentrations n k (t) within the numerical precision.
The structure of the kinetic kernels (7) allows to apply highly efficient numerical methods, in particular, the fast and accurate method of time-integration of Smoluchowski-type equations [35][36][37][38][39]. The efficiency and accuracy of this approach in solving the aggregating-andshattering equations has been demonstrated in Ref. [30], where the numerical results have been compared with the available analytical solutions [10].

A. Steady state size distribution
Solving numerically Eqs. (3)-(4) with kernel (7) for a < 1/2, we observe that the concentrations relax monotonically to a steady-state, see Fig. 2. In Fig. 2 we also compare the numerical results with the analytical solution for the steady-state distribution n k , discussed below. The numerical and analytical solutions agree fairly well.
Similar behavior is observed for the general kernel (5). When θ = µ − ν < 1, the concentrations relax monotonically to a steady-state, and the final distribution agrees with the one predicted theoretically, see Fig. 3. The steady-state size distribution may be interpreted as a stable fixed point in the language of dynamical systems [40].
For a ≥ 1/2 a relaxation to a steady-state distribution occurs through oscillations, provided the parameter λ, quantifying the shattering intensity, is relatively small. This is illustrated in Fig. 4, where the time dependence of the total number of aggregates, N (t) = k≥1 n k (t), is shown; the figure also demonstrates that the oscillations are more pronounced and persist for longer time as a increases, while λ decreases.
We found the oscillations independently of initial conditions; here we use the mono-disperse initial conditions, n k (0) = M δ 1,k and step-wise initial conditions with the same total mass M = 5.5. Unless explicitly stated, the reported results refer to the initial conditions (9). For a → 1 and relatively small λ we observe stable, seemingly never-ending oscillations, see Fig. 5, where the temporal behaviors of the total density N (t) and the second moment M 2 (t) = k≥1 k 2 n k (t) are depicted.
Making the time averaging of the densities over the oscillation period, one obtains the distribution of the averaged quantities n k osc , which has a form of the powerlaw with a cutoff at k ∼ k 0 , see Fig. 6: 2. General kernels (5) We observed oscillations for the general kernel (5) when θ = ν − µ > 1 (which corresponds to a > 1/2 of the Brownian kernel). However, if θ is not close to θ = 2, the system relaxes to a steady distribution through the damped oscillations, even for rather small λ, see Fig. 7.
Our simulations imply the existence of a critical value λ c (θ) such that for λ < λ c (θ) in the long time limit the  (7) with a = 0.9 and different λ. Seemingly neverending oscillations are observed for λ = 0.005.
system approaches a limit cycle, viz. concentrations exhibit never-ending oscillations. This has been checked for the Brownian kernel and for the general kernel (5). Although for a < 0.9 we have observed only damped oscillations, we believe that never-ending oscillations would emerge for all a > 1/2 and sufficiently small λ. This is seemingly true for the general case: The steady oscillations would be observed for any θ > 1 if λ is small enough. We cannot prove this numerically due to unaccessible number of equations needed to simulate the systems with such small λ. For instance, to simulate the system with a = 0.9 and λ = 0.005 depicted in Fig. 5, more than 250 000 equations have been used. Our estimates (discussed in Ref. [30]) indicate that the number of equations N eq , needed to guarantee a requested accuracy rapidly grows with the decreasing λ. To simulate a system with λ < λ c for a < 0.9 one needs more than a millions equations which is too large for practical implementation. Nevertheless, based on our results, we formulate the following Conjecture. (i) When θ = µ − ν < 1, the system has a single stable fixed point for all values of λ; the steady state distribution of cluster sizes n k corresponds to this stable point. (ii) When θ > 1, there exists a critical λ c (θ), such that for λ ≥ λ c the system possesses a stable fixed point with the according distribution n k . This may be a stable focus for some values of λ manifesting in damped oscillations. (iii) When θ > 1 and λ < λ c (θ), the system possesses a stable limit cycle.
As it follows from our numerical results, the critical shattering λ c strongly depends on the exponent θ; its dependence on the other exponent β = µ + ν seems to be weak (if any), but is still to be studied.

III. THEORETICAL ANALYSIS
To explain theoretically the observed behavior of the aggregation-and-shattering systems we analyze separately the systems that attain a steady-state distribu-tion and those that demonstrate never-ending oscillations. For the former case we apply the asymptotic analysis, while in the latter situation we analyze oscillations qualitatively.
A. Asymptotic analysis of a steady-state cluster size distribution In Ref. [30] we gave a condensed account of the derivation of the steady-state distribution; here we present a more detailed derivation.
When the system reaches a steady-state, differential equations (3)-(4) become algebraic equations To analyze these equations we introduce the generating functions and the moments Multiplying (11) by z k and summing over all k ≥ 1 we arrive at Using C γ (1) = M γ and specializing (13) to z = 1 we obtain To analyze n k for k 1 we will use the above equations and exploit standard methods of asymptotic analysis to extract the behavior of the generation functions C γ (z). We consider separately kernels with θ < 1 and θ > 1.
Let us assume that for k 1 our steady-state distribution has a similar form: with yet unknown τ , ω and C. Equation (15) implies where z 0 = e ω and z = z/z 0 . Obviously, C γ (z ) diverges for z > 1 and converges for z < 1 for all γ and τ . We assume, that C γ (z = 1) exists, that is, k≥1 k γ−τ converges. The tail of n k is reflected in the behavior of C γ (z ) when z → 1 − 0. Suppose k≥1 k γ−τ +1 diverges. Still, k≥1 k γ−τ +1 (z ) k converges for z < 1. The closer z is to 1, the larger the size of the clusters k, that make the main contribution to C γ (z ). Hence the dependence of C γ (z ) on z for z → 1 characterizes the dependence of n k on k for k 1. To quantify this relation we differentiate C γ (z) with respect to z: is the gamma-function and we use z → 1 −0.
Integrating with respect to z we obtain Substituting C γ (z) with γ = ν and γ = µ into Eq. (13) we obtain terms with different powers of (1 − z ). To satisfy this equation we equate to zero all these terms separately. The zero-order terms yield The terms of the order (1 − z ) τ +γ−1 with γ = ν and γ = µ imply Finally, the rest of the terms should satisfy from which 2τ − ν − µ − 2 = 1, or τ = 3 + β 2 (22) where β = ν + µ. Now we substitute which follows from (19) and (20) into (18) to obtain From Eqs. (24) and (14) we get We have ω λ 2 − 2λ 3 + . . . λ 2 for small λ leading to To estimate the constant C we utilize the distribution (26) together with mass conservation to yield In the non-gelling β < 1 region, the major contribution to the integral in (27) comes from k 1, so the usage of (26) is justified. For the Brownian kernel ν = −µ = a and β = 0, so the amplitude is C λM/ √ π and n k λM √ πk 3/2 e −λ 2 k for k 1.

B. Qualitative analysis
To understand the mechanism of the stable oscillations let us consider a special Brownian kernel with ν = −µ = a = 1. The monomer density satisfies (30) and the rate equation for the cluster density is (We set M = 1). These equations are not closed as they involve the moments M −1 and M 2 . One can write the rate equation for M 2 , but it involves the third moment M 3 . This continues ad infinitum leading to analytically unsolvable hierarchy. Note that on the right-hand side of Eq. (30) one term is negative and the other one is positive; the first is of the order of n 1 and the second is of the order λM 2 . Initially only small clusters present in the system, so that λM 2 is small for λ 1 and n 1 decreases. Due to the conservation of mass the decrease of n 1 implies the increase of other concentrations. Hence, after some time a wider cluster size distribution is established, such that λM 2 increases. When it exceeds n 1 , the right-hand side of Eq. (30) becomes positive and n 1 starts to grow. Due to the conservation of mass, the growth of n 1 implies the decay of other concentrations, which leads to the decrease of M 2 and eventually to the negative sign of the right-hand side of Eq. (30). Then the cycle repeats.
Let us try to put the above narrative picture into somewhat more quantitative terms. Firstly, we notice that the oscillations of concentrations correspond to the periodically varying distribution of cluster sizes, as it is illustrated in Fig. 9. Roughly speaking, the size distribution n k (t) behaves in such a way that the effective slope of this distribution α(t) and the effective cutoff k max (t) periodically change in time. Averaging over these oscillations we obtain the distribution n k osc depicted in Fig. 6.
To understand the nature of the observed behavior of the system, we develop a qualitative theory. To this end we approximate the real distribution n k (t) by a model distribution n mod k (t) which reflects the most prominent features of the real distribution. Namely we assume that n k (t) may be characterized by a power-law distribution with a varying slope α(t) and a large-size cutoff k max (t); it should obey the mass conservation. Namely, we assume the following model distribution: The real distribution of the aggregate sizes n k (t) may be approximated by the model distribution (32) applying a coarse graining. For the qualitative analysis addressed here we exploited the most simple approach. Namely, we use the value of n 1 (t), obtained in the simulations, and find the parameters α(t) and k max (t) from the conservation of mass. That is, we numerically find the pair (α, k max ) with integer k max , that minimizes the difference, the generalized harmonic numbers. The variation with time of the model distribution (32), which mimics the real distribution, is shown in Fig. 9. In Fig. 10 we show the periodic variation on the model parameters α(t) and n 1 (t) and demonstrate that the variation of the slope is limited by the interval 1 < α(t) < 2.
To perform a qualitative analysis we focus on the qualitative dependence of the moments M −1 , N , M 2 on n 1 , α, k max and approximate the summation by integration: where we introduce the coefficients b i . These coefficients (assumed to be constant) are of the order of one and account for the difference between integration and summation. Hence we obtain where we have used the condition k max 1. Similarly, the conservation of mass yields the relation between n 1 (t), α(t), and k max (t): Using (34)-(35), we recast (30)-(31) into To show that never-ending oscillations are possible we perform the linear stability analysis of Eqs. (36a)-(36b). We consider the coefficients b i (i = −1, 0, 1, 2) as known and of the order of unity. Further, we assume that there is a fixed point, n 1 = n (0) 1 and α = α 0 . At the fixed point G 1 (n n 1 , α) and for G w (n 1 , α) denote the right-hand sides of (36a) and (36b), respectively. We also shortly write and deduce the linearized equations d dt for the deviations δn 1 = n 1 − n (0) 1 and δα = α − α 0 . The eigenvalues of the matrix in (38) are Oscillations may occur if the above eigenvalues possess an imaginary part. This condition, Im(ν 1/2 ) = 0, require the negatives determinant in If the real part of the eigenvalues is negative, that is (g 1n + g 2α )/2 < 0, the fixed point is stable; in this case the cluster distribution relaxes to the stead-state n k . In the opposite case of positive real part, (g 1n + g 2α )/2 > 0, the fixed point is linearly unstable and the oscillations grow and eventually stabilized by non-linear terms.
The coefficients b i are unknown, so we cannot locate the fixed point (n 1 ≤ 10λ and 1 < α 0 < 1.5 for different λ. The results are shown in Fig. 11. Figure 11 demonstrates that for λ = 0.0001 there is a large area in the domain of the α 0 , n (0) 1 plane where steady oscillations may be observed. These may be either linearly stable oscillations or the growing ones, stabilized by non-linear terms. For relatively large λ = 0.1, the steady oscillations may arise only in a tiny part of the α 0 , n (0) 1 plane. This corresponds to the kinetic regimes observed for the full set of aggregation-fragmentation equations: the emergence of the oscillations for small values of λ and their absence for large λ.
We also note that the average slope of the concentration distribution, α 0 , is located within the interval 1.05 ≤ α 0 ≤ 1.45, see Fig. 11, with the median of 1.25. This is consistent with the slope of the averaged over oscillations distribution n k osc depicted in Fig. 6.

C. Concentration oscillations in thermodynamically closed systems
As noted above, the time-independent rates K ij and F ij = λK ij imply a steady supply of energy. This follows generally from the second law of thermodynamics, which excludes steady cyclic processes without energy supply and may be illustrated for a particular microscopic mechanism of a ballistic aggregation and shattering; this happens in planetary rings or atmospheric processes. Indeed, the conservation of momentum of coalescing particles dic-tates a withdrawal of a part of their kinetic energy, associated with the relative motion. This energy is transmitted either to the internal degrees of freedom of the particles (as in planetary rings), or to the surrounding gas (as in atmospheric processes). Similarly, the total kinetic energy of fragments is smaller than the initial kinetic energy of the colliding aggregates, since part of the energy is spent to break inter-fragment bonds. Hence, both aggregation and fragmentation processes lead to a gradual reduction of the total kinetic energy of the system. This causes a slowdown of the both processes, and the respective decrease of the rate coefficients. Here we consider thermodynamically closed systems, where the energy supply is lacking. Namely, we consider systems of particles undergoing ballistic aggregation and fragmentation. We chose such systems since the corresponding microscopic rates K ij and F ij are available [10,13,16], see the Appendix.
Physically, the decay of kinetic energy causes a permanent decrease of collision frequency and the decrease of the kinetic rates. Moreover, the fragmentation rates additionally decrease, since the fraction of fast particles, which cause shattering, also drops down. Referring for detail to the Appendix, we present here the equations for aggregation and shattering for thermodynamically closed systems: while for heavier clusters, k ≥ 2, we have K ij appearing in Eqs. (41)- (42) are defined by (5). Time is measured in collision units, τ −1 c = 2 √ 2πσ 2 n 1,0 T (t)/m, where σ, m and n 1,0 are respectively the diameter, mass and initial concentration of monomers and T (t) is the characteristic temperature. The concentration of aggregates is measured in units of the initial concentration of monomers, n 1,0 . The constants A, B and b in Eq. (43) are expressed respectively in terms of the characteristic fragmentation energy and temperature decay rate (see Appendix for the detail).
The results of numerical solution of Eqs. (41)-(43) are presented in Fig. 12. As may be seen from the figure the persistent concentration oscillations emerge in thermodynamically closed systems. Since the time for the depicted oscillations is measured in the collision units, one concludes that the period of these oscillations steadily increases with time; this is also visible in the collision time scale. As one can see from Fig. 12, in thermodynamically closed systems there exist a regime, when the oscillations first decay and then again grow.

IV. CONCLUSIONS
We have studied numerically and analytically a class of aggregation-fragmentation models. Mathematically, the problem is described by an infinite set of Smoluchowskilike equations with the homogeneous aggregation and fragmentation kernels which respectively read K i,j = i ν j µ + j ν i µ and F i,j = λK i,j , where the parameter λ quantifies the intensity of fragmentation. We consider the case of a complete decomposition (shattering) of colliding aggregates into monomers. This model and a similar model, with a source of monomers and evaporation (instead of shattering) of large clusters has been studied recently in [17,18]. For the kernels with θ = ν − µ < 1 we obtain an analytical solution for the steady-state size distribution of the aggregates n k and confirm numerically the relaxation of the size distribution to this steady-state form. For kernels with θ = ν − µ > 1, we observe that the dynamic of the system dramatically depends on the value of the fragmentation constant λ. While for λ < λ c the system relaxes to a steady-state through damped oscillations of concentrations, for λ ≥ λ c no steady-state distribution of the cluster size has been detected.
The emergence of stable oscillations in a closed system of aggregating and fragmenting particles, that lacks any sinks and sources of mass, and formally corresponds to an infinite number of species, is new and surprising. Persistent oscillations have been detected not only for systems, closed with respect to the total mass, but also for thermodynamically closed systems, when the notion "thermodynamics" is meaningful. In Ref. [27] stable oscillations have been detected numerically for Smoluchowski equations for an open system of reversibly aggregating particles (without fragmentation) with a source of monomers and sink of large clusters, which makes the system finite. For a small closed system comprising monomers, dimers, trimers and exited monomers, stable oscillations of concentrations have been also reported [28]. Similarly, steady chemical oscillations have been found in a simple dimerization model (see e.g. [29] and references therein).
Our findings may help to understand various phenomena observed in the systems with aggregation and fragmentation, in particular the periodic formation and destruction of clumps in F Ring of Saturn [41], where particles of different mass suffer aggregative and disruptive impacts, presumably under the mass conservation condition. A complete understanding of this phenomenon is presently lacking.

Appendix A
General expressions for the aggregation and fragmentation rates for a system of ballistically moving particles (molecules, macroscopic grains, etc.) that suffer pairwise collisions, have been reported in [10,13]. These rates read in present notations. Here T i are partial temperatures of aggregates of size i, and mass m i = m 1 i (m 1 is the mass of monomer), that characterizes the average kinetic energy of such aggregates [10,13,16]. E agg and E frag are respectively the aggregation and fragmentation energy,Ẽ agg = E agg /ε 2 , where ε is the coefficient of normal restitution that characterizes dissipative losses in the impacts [16]. σ ij = σ 1 i 1/D + j 1/D , where σ 1 is monomer diameter and D is the dimension of the aggregates, which may be fractal. Based on the results of Ref. [42], we assume that partial temperatures scale as T i (t) = T (t)i γ , where T (t) is the characteristic temperature of the gas mixture. (It may be shown that after a short relaxation time the rate of change of temperatures of all species is the same, T −1 i dT i /dt = T −1 dT /dt, [42]). It is convenient to recast the above kinetic coefficients into the fform where n 1,0 is the initial concentration of monomers, which we will use as a unit of concentration, T 0 is the initial characteristic temperature, and gives the initial characteristic collision frequency. The quantity appearing in (A2) is the effective average ratio of fragmentation and kinetic energy. We also assume for simplicity that the aggregation energy is large, so that B ijẼagg 1. The dimensionless kernelsK ij are obtained by a straightforward solution of the Boltzmann equation [10,13,16]. Here we apply a standard simplification [1,2] for these kernels, which allows analytical analysis. It is based on the observation that the main properties of the solutions to the Smoluchowski equations depend on two indices β and β 1 , characterizing the kinetic rates K ij . The first index quantifies the homogeneity degree of a kernel, and the second one the size dependence at the maximal size asymmetry. Namely, K ai,aj ∼ a βK ij ;K 1,j ∼ j β1 for i, j 1.
For the kernelsK ij = i ν j µ + i µ j ν introduced in Eq. (5), one obtains β = ν + µ and β 1 = max(ν, µ). Therefore Eqs. (A2) show that ν and µ are related to the physical parameters D and γ via relations (ν, µ) = (2/D, (γ − 1)/2) for γ < 1 and (ν, µ) = (2/D + (γ − 1)/2, 0) for γ > 1. Next, we derive the equation for the characteristic temperature T (t). This may be done using the approach of Ref. [16], which yields Here Q ij (T ) are temperature-dependent rate coefficients and Γ i describes the energy input to the system due to the interaction of the aggregates of size i with the external sources of energy (see also [42]). Here we do not need explicit expressions for these quantities. We just state that the presence of the energy sources Γ i in Eq. (A4) yields the solutions with a constant temperature T = const., corresponding to the systems with time-independent rates K ij and F ij . For thermodynamically closed systems, the temperature commonly decreases with time (see the discussion in Ref. [16]). The solutions of the coupled Smoluchowski-like equations (A2) and (A4) for concentrations and temperature is beyond the scope of the present study. For the qualitative analysis, we assume a power-law decay of the characteristic temperature with time, T = T 0 (1 + t/τ 0 ) −δ ; such assumption is justified by the results of Ref. [16]. The value of δ depends on the parameters of the system and may vary in a wide interval [16]. Using the collision frequency, τ −1 c (t) = τ −1 0 [T (t)/T 0 ] 1/2 at the current time t, we introduce a new dimensionless timet, measured in collision units. It is related to the laboratory time as τ −1 c (t)dt = dt. The dependence of temperature on the new time then reads where b = 2δ/(2 − δ) and B is a constant. The kinetic rates may be expressed in terms of the collision-based timet: In our simulations we choose b = 0.2 (which corresponds to δ = 10/11). Substituting the rates K ij and λ from Eqs. (A6) and (A7) into Eqs. (3) and (4)