Oscillations in aggregation-shattering processes
S. A. Matveev, P. L. Krapivsky, A. P. Smirnov, E. E. Tyrtyshnikov, N. V. Brilliantov
OOscillations in aggregation-shattering processes
S. A. Matveev , , , P. L. Krapivsky , A. P. Smirnov , , E. E. Tyrtyshnikov , , and N. V. Brilliantov Skolkovo Institute of Science and Technology, Moscow, Russia Faculty of Computational Mathematics and Cybernetics, Lomonosov MSU, Moscow Russia Institute of Numerical Mathematics RAS, Moscow Russia Department of Physics, Boston University, Boston, MA 02215, USA and Department of Mathematics, University of Leicester, Leicester LE1 7RH, United Kingdom (Dated: September 17, 2018)We observe never-ending oscillations in systems undergoing aggregation and collision-controlledshattering. Specifically, we investigate aggregation-shattering processes with aggregation kernels K i,j = ( i/j ) a + ( j/i ) a and shattering kernels F i,j = λK i,j , where i and j are cluster sizes andparameter λ quantifies the strength of shattering. When 0 ≤ a < /
2, there are no oscillationsand the system monotonically approaches to a steady state for all values of λ ; in this region weobtain an analytical solution for the stationary cluster size distribution. Numerical solutions ofthe rate equations show that oscillations emerge in the 1 / < a ≤ λ is sufficientlylarge oscillations decay and eventually disappear, while for λ < λ c ( a ) oscillations apparently persistforever. Thus never-ending oscillations can arise in closed aggregation-shattering processes withoutsinks and sources of particles. Two complimentary processes, aggregation and frag-mentation [1–3], occur in numerous systems that dra-matically differ in their spatial and temporal scales. Re-versible polymerization in solutions [1] and merging ofprions (cell proteins) [4] are typical examples on themolecular scale. On somewhat larger scales airborneparticles perform Brownian motion in atmosphere andcoalesce giving rise to smog droplets [5]. Aggregationof users in the Internet leads to the emergence of com-munities and forums [2, 6] which can further merge orsplit. Vortexes in a fluid flow merge and decomposeforming turbulent cascades [7]. On much larger scales,aggregation-fragmentation processes take place in plane-tary rings, like Saturn rings, where the particle size dis-tribution is determined by a subtle balance between ag-gregation and fragmentation of the rings particles [8–12].In spatially homogeneous well-mixed systems, aggre-gation and fragmentation processes are described byan infinite set of nonlinear ordinary differential equa-tions (ODEs) for the concentrations of clusters of var-ious masses. Such equations are intractable apart froma few special cases. The long-time behavior, however,is occasionally known—the processes of aggregation andfragmentation act in the opposite directions and hencethe cluster size distribution often becomes stationary inthe long time limit [34]. The emergence of the station-ary cluster size distribution can be mathematically inter-preted as the manifestation of the fixed point of the dif-ferential equations [13]. For a single differential equation,fixed points determine the long time behavior, while fortwo coupled differential equations the asymptotic behav-ior may be determined by a fixed point or a limit cycle. Inthe case of infinitely many coupled ODEs, limit cycles arefeasible, yet they haven’t been observed in aggregation-fragmentation processes. More precisely, there were signsof oscillations in a few open systems usually driven byconstant source of monomers and by sink of large clus-ters. In this paper we report oscillations in a closed sys- tem undergoing aggregation and collision-controlled frag-mentation.In the most important case of binary aggregationthe collision between two clusters comprising i and j monomers may result in the formation of a joint aggre-gate of i + j monomers. Symbolically, [ i ]+[ j ] K i,j −−−→ [ i + j ] , where K ij is the merging rate (see Fig. 1). Let n k bethe concentration of clusters that contain k monomers.These quantities obey the Smoluchowski equations [2, 3]: dn k dt = 12 (cid:88) i + j = k K i,j n i n j − n k ∞ (cid:88) i =1 K i,k n i . (1)The first gain term on the right-hand side gives the for-mation rate of k -mers from smaller clusters, while thesecond terms describes the disappearance of k -mers dueto collisions with other clusters (the factor 1 / ji i ij K j a i ji ij F j b FIG. 1: Aggregation (a) and shattering (b) of clusters.
In this article we consider collision-controlled fragmen-tation, which is thought to be responsible e.g. for inter-stellar dust clouds and planetary rings [8, 9, 11, 14]. Weexplore the extreme version, namely a complete shatter-ing of two colliding partners into monomers (see Fig. 1).Symbolically [ i ] + [ j ] F i,j −−→ [1] + [1] + . . . [1] (cid:124) (cid:123)(cid:122) (cid:125) i + j where F i,j quantifies the shattering rate. It has been shown [8] thatthis shattering model is rather generic—realistic impact a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug models with a strong dominance of small debris over thelarge ones yield the same resulting cluster size distribu-tion n k . Following [8], we assume that the shattering andaggregation kernels are proportional, F i,j = λK i,j . (2)The parameter λ characterizes the relative frequency ofcollisions leading to merging and shattering.Incorporating the shattering process with the shatter-ing kernel (2) into Eqs. (1) we arrive at dn dt = − n ∞ (cid:88) i =1 K ,i n i + λ ∞ (cid:88) i =2 ∞ (cid:88) j =2 ( i + j ) K i,j n i n j + λn ∞ (cid:88) j =2 jK ,j n j (3) dn k dt = 12 k − (cid:88) i =1 K i,k − i n i n k − i − (1 + λ ) n k ∞ (cid:88) i =1 K k,i n i . Shattering leads to the increase of monomers explainingthe gain terms in the first equation (3) and it leads tothe decrease of the density of other clusters explainingthe loss term in the second equation.A microscopic analysis is needed to establish how thekernels K i,j and F i,j depend on the masses, see e.g. [8,11]. The kernels are always symmetric, K i,j = K j,i , andin most applications homogeneous functions of i and j .Aggregation-shattering equations (3) for the generalizedproduct kernels, K i,j = ( ij ) µ , have been investigated in[8]. A more general family of kernels, K i,j = i ν j µ + i µ j ν ,is often used in studies of aggregation, see e.g. Ref. [15]where a source of monomers and sink of large clusterswas present. We shall focus on a special case of µ = − ν , K i,j = i a j − a + i − a j a , (4)which is also known as a generalized Brownian kernel [16].Below we always assume that a ≤
1, since aggregationequations with kernel (4) satisfying a > a = 0). The steady-state solutions have beenobtained for a wider class of models, including irre-versible aggregation model with a monomer source [23],aggregation-fragmentation model with the generalizedproduct kernel [8] and for a somewhat similar open sys-tem with a source of monomers and collisional evapora-tion of clusters with the kernel K i,j = i ν j µ + i µ j ν [15].An open aggregating system with the same coagulationkernel driven by input of monomers and supplementedwith removal of large clusters has been studied in [24].Stable oscillations have been numerically observed [24]in this system with finite number of cluster species. For a closed system consisting of monomers, dimers, trimersand exited monomers, stable oscillations of concentra-tions have been reported [25]. Steady chemical oscilla-tions have been also found in a simple dimerization model(see e.g. [26] and references therein).Here we consider closed systems undergoing aggrega-tion and shattering processes, so there are no sourcesand sinks of monomers and clusters. The aggregationand shattering kernels are described by Eqs. (4) and (2).One expects that in the closed system with two oppositeprocesses and without sinks and sources, a steady stateis achieved. This is indeed the case when the parameter a < /
2. Surprisingly, for 1 / < a ≤ λ , a steady state is not reached and instead clusterconcentrations undergo never-ending oscillations.An important advantage of the kernel (4) is the possi-bility to apply highly efficient numerical methods. Herewe exploit a fast and accurate method of time-integrationof Smoluchowski-type kinetic equations developed in re-cent studies [27–32], that has been adopted for the dis-crete distribution of cluster sizes. In our simulations weuse up to N eq = 2 ≡ ,
288 equations; in practice,we choose N eq in such a way, that the further increase of N eq does not impact the simulation results for n k withinthe numerical accuracy [35].We confirm the efficiency and accuracy of the abovefast-integration method for constant kernels ( a = 0),comparing the numerical results with the available an-alytical solutions [8] and find that the smaller the pa-rameter λ the longer it takes for the system to reach thesteady state, see the Supplementary Material (SM). Thistendency persists for the kernels (4) with a >
0. Wealso observe that for a < / n k ( t ). Moreover, the steady-state distributionof the cluster concentrations agrees fairly well with ana-lytical results for n k derived below. Figure 2 illustratesthe numerical solution for steady-state distribution for a = 0 .
05 and a = 0 . λ . In the language of dynamical systems[13] we conclude that the system possesses a stable fixedpoint with the steady-state cluster size distribution.For 1 / < a ≤
1, we also observed the relaxation to asteady state for sufficiently large λ , the relaxation occursthrough oscillations for smaller λ , and when λ < λ c ( a )the oscillations persist. We detected oscillations indepen-dently of the initial conditions.The dynamic of the system described by Eqs. (3) and(4) is invariant with respect to re-scaling of the total massdensity M = (cid:80) ∞ k =1 kn k (see the SM). Below we reportsimulation results for the stepwise initial distribution n k (0) = (cid:40) . k = 1 , . . . , k >
10 (5)and we also simulated the evolution starting with mono-disperse initial condition, n k (0) = M δ ,k , with the same -12 -8 -4 n k kAnalyticalNumerical -12 -8 -4 n k kAnalyticalNumerical FIG. 2: Comparison of the steady-state numerical solution ofEqs. (3) for the kernel (4) with a = 0 . λ = 0 .
003 (top) and a = 0 . λ = 0 .
02 (bottom) with the analytical steady-statesolution, Eq. (9). mass M = 5 .
5. Unless explicitly stated, the results belowcorrespond to the initial condition (5).In Figs. 3–4 we demonstrate the time dependence ofthe cluster density N ( t ) = (cid:80) ∞ k =1 n k ( t ). Figure 3 showsthat in the range 0 . ≤ a ≤ . . ≤ λ ≤ .
01, theoscillations become more pronounced when a increasesand λ decreases.In Fig. 4 we show oscillating solutions for N ( t ) for0 . < a ≤
1. The new feature observed in Fig. 4 is theemergence of stationary oscillations. All cluster concen-tration n k ( t ) perform these stable oscillations; the formof the oscillations depends on the cluster size and theamplitude decreases with the increasing size, see Fig. 5.Figure 6, demonstrates that the system reaches a limitcycle [36] which does not depend on the total mass orinitial conditions.Our results indicate the existence of a critical value of λ c ( a ) such that for λ < λ c ( a ) the steady-state solutionis no longer stable and instead the system approaches toa limit cycle. Although for a < . a ≥ / λ are too small. When λ is small, a huge numberof equations is needed to achieve a requested precision.For example, for the simulations presented in Fig. 3 al-ready 200 ,
000 equations have been used. For systemswith λ < λ c for a < .
9, one needs to solve N eq > nonlinear ODEs which is a formidable task even whenfast numerical methods are applied (see the SM).Our major observations may be summarized as follows: N ( t ) time λ = 0.005 λ = 0.003 λ = 0.001 N ( t ) time λ = 0.01 λ = 0.002 λ = 0.001 FIG. 3: Time dependence of the clusters density, N ( t ), for a = 0 . a = 0 .
75 (bottom) and different λ . For allthese systems damped oscillations are found that tend to asteady-state. The oscillations become more pronounced withincreasing a and decreasing λ .
1. When a < /
2, there exists a single stable fixedpoint for all values of λ ; the steady state distri-bution of cluster sizes n k corresponds to this fixedpoint.2. When 1 / < a ≤ λ > λ c ( a ), the system hasa single stable fixed point; it may be a stable focus,resulting in dumped oscillations.3. When 1 / < a ≤ λ < λ c ( a ), the system hasan attractive limit cycle.In the 1 / < a ≤ dn /dt = 0 and dn k /dt = 0 in Eqs. (3) and solve the re-sulting infinite system of algebraic equations. Introduc-ing the generating functions C ± a ( z ) = (cid:80) k ≥ k ± a n k z k ,we transform (3) into C a ( z ) C − a ( z ) + (1 + λ ) zn ( M a + M − a ) (6)= (1 + λ ) ( M a C − a ( z ) − M − a C a ( z )) . where C ± a (1) = M ± a . Setting z = 1 in Eq. (6) yields M a M − a = 1 + λ λ n ( M a + M − a ) . (7)To analyze the tail of the size distribution, i.e. n k for k (cid:29)
1, we exploit the methods described e.g. in [2, 16] N ( t ) time λ = 0.005 λ = 0.01 λ = 0.02 FIG. 4: Time dependence of the cluster density, N ( t ) for a = 0 . λ . For small λ < λ c ( a ) stable oscillationsemerge. -12 -8 -4
0 50 100 150 200 n k timen n n n -12 -8 -4
0 50 100 150 200 n k time FIG. 5: Stationary oscillations of the aggregate concentra-tions for a = 0 .
95 (top) and a = 1 (bottom) and λ < λ c ( a ).The shape of the oscillations depends on the aggregates’ size;the amplitude of the oscillations decreases with the size. in the context of similar problems. Recalling that when a = 0 the tail is n k (cid:39) λπ − / k − / e − λ k for k (cid:29) n k (cid:39) Ck − τ e − ωk for k (cid:29)
1, withsome constants C , τ and ω . Expanding the generatingfunctions C ± a ( z ) near z (cid:48) → −
0, where z (cid:48) = z/z and z = e ω , we get C ± a ( z ) = C ± a ( z ) + C Γ(1 ± a − τ )(1 − z (cid:48) ) τ ∓ a − . (8)Here Γ( x ) is the gamma function and we assume that C a ( z ) = (cid:80) k ≥ k a n k z k < ∞ exists for given a and τ . Ifwe substitute the above C ± a ( z ) into Eq. (6) we obtainterms with different powers of the factor (1 − z (cid:48) ). Tosatisfy this equation we equate to zero all these terms n n Monodisperse,M=3Polydisperse,M=3Monodisperse,M=5.5Polydisperse,M=5.5
FIG. 6: Limit cycle for the steady-state oscillations in termsof n ( t ) and n ( t ) for a = 0 .
95 and λ = 0 . M = 5 . n k = Mδ ,k and stepwise, (5). For M = 3 the initialand current values of n ( t ) and n ( t ) are re-scaled accordingly.The relaxation to the unique limit cycle is clearly visible. separately. In this way we obtain equations for the zero-order terms of (1 − z (cid:48) ), and the terms of the order of(1 − z (cid:48) ) τ ± a − . Combining these equations with Eq. (7)we find τ = 3 / ω (cid:39) λ , and finally the amplitude C = M λπ − / (see SM for details). Thus the tail of thecluster size distribution reads n k (cid:39) λπ − / M k − / e − λ k for k (cid:29) . (9)In the above analysis we assume that C a ( z ) exists. Thisis a consistent assumption when a < /
2, but fails for a ≥ / K i,j = ( i/j ) a + ( j/i ) a and shattering kernel F i,j = λK i,j . For a < /
2, we obtained an analyticalsolution for the steady-state cluster size distribution andconfirmed numerically the relaxation of the size distribu-tion to this steady-state form. For a ≥ /
2, the temporalbehavior drastically depends on the shattering constant λ : When λ > λ c ( a ) the system relaxes to a steady-statethrough damped oscillations, while for λ < λ c ( a ) oscilla-tions become stationary and persist forever (Figs. 4–6).Using the language of dynamical systems our observa-tions can be reformulated as follows: (i) for a < / λ , (ii) for 1 / ≤ a ≤
1, the systemhas a single stable fixed point (which may be a stablefocus) when λ ≥ λ c ( a ), and (iii) for 1 / ≤ a ≤ λ < λ c ( a ) the system possesses a stable limit cycle.Limit cycles may arise already for two coupled ODEs[13]. Still, the emergence of stable oscillations in a closed system comprising an infinite number of species and un-dergoing aggregating and shattering is striking. To thebest of our knowledge this phenomenon has not beenpreviously observed and a relaxation towards the steady state was believed to be the only possible scenario.The work was supported by the Russian Science Foun-dation, grant 14-11-00806. [1] P. J. Flory, Principles of Polymer Chemistry (CornellUniversity Press, 1953).[2] P. L. Krapivsky, S. Redner, and E. Ben-Naim,
A kineticview of statistical physics (Cambridge University Press,2010).[3] F. Leyvraz, Physics Reports , 95 (2003).[4] T. Poeschel, N. V. Brilliantov, and C. Frommel, Biophys.J. , 3460 (2003).[5] R. C. Shrivastava, J. Atom. Sci. , 1317 (1982).[6] S. N. Dorogovtsev and J. F. F. Mendes, Evolution ofnetworks: From biological nets to the Internet and WWW (Oxford University Press, 2003).[7] V. E. Zakharov, V. S. L’vov, and G. Falkovich,
Kol-mogorov Spectra of Turbulence I: Wave Turbulence (Springer, 2012).[8] N. V. Brilliantov, P. L. Krapivsky, A. Bodrova, F. Spahn,H. Hayakawa, V. Stadnichuk, and J. Schmidt, PNAS , 9536 (2015).[9] V. Stadnichuk, A. Bodrova, and N. V. Brilliantov, Int. J.Mod. Phys. B , 1550208 (2015).[10] J. N. Cuzzi, J. A. Burns, S. Charnoz, R. N. Clark, J. E.Colwell, L. Dones, L. W. Esposito, G. Filacchione, R. G.French, M. M. Hedman, et al., Science , 1470 (2010).[11] N. V. Brilliantov, A. Bodrova, and P. L. Krapivsky,J. Stat. Mech.: Theory and Experiment , P06011(2009).[12] L. Esposito, Planetary Rings (Cambridge UniversityPress, 2006).[13] S. H. Strogatz,
Nonlinear Dynamics and Chaos (AddisonWesley publishing company, 1994).[14] P. L. Krapivsky and E. Ben-Naim, Phys. Rev. E ,021102 (2003).[15] C. Connaughton, A. Dutta, R. Rajesh, and O. Zaboron-ski, EPL , 10002 (2017).[16] P. L. Krapivsky and C. Connaughton, J. Chem. Phys. , 204901 (2012).[17] E. M. Hendriks, M. H. Ernst, and R. M. Ziff, J. Stat. Phys. , 519 (1983).[18] P. G. J. van Dongen, J. Phys. A , 1889 (1987).[19] N. V. Brilliantov and P. L. Krapivsky, J. Phys. A: Math.Gen. , 4787 (1991).[20] P. Lauren¸cot, Nonlinearity , 229 (1999).[21] L. Malyushkin and J. Goodman, Icarus , 314 (2001).[22] R. C. Ball, C. Connaughton, T. H. M. Stein, andO. Zaboronski, Phys. Rev. E , 011111 (2011).[23] H. Hayakawa, J. of Phys. A , L801 (1987).[24] R. C. Ball, C. Connaughton, P. P. Jones, R. Rajesh, andO. Zaboronski, Phys. Rev. Lett. , 168304 (2012).[25] V. I. Bykov and A. N. Gorban, Chem. Eng. Sci. , 1249(1987).[26] M. Stich, C. Blanco, and D. Hochberg, Physical Chem-istry Chemical Physics , 255 (2013).[27] S. A. Matveev, A. P. Smirnov, and E. E. Tyrtyshnikov,J. Comput. Phys. , 23 (2015).[28] A. P. Smirnov, S. Matveev, D. Zheltkov, E. E. Tyrtysh-nikov, et al., Procedia Computer Science , 2141 (2016).[29] S. A. Matveev, D. A. Zheltkov, E. E. Tyrtyshnikov, andA. P. Smirnov, J. Comput. Phys. , 164 (2016).[30] A. Chaudhury, I. Oseledets, and R. Ramachandran,Computers & Chemical Engineering , 234 (2014).[31] W. Hackbusch, Computing , 145 (2006).[32] W. Hackbusch, Numerische Mathematik , 627 (2007).[33] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E ,061132 (2008).[34] There are a few exceptions when the typical cluster sizediverges and/or the system undergoes a non-equilibriumphase transition, see e.g. [33].[35] In the Supplementary Material (SM) we justify that aninfinite system (3) with the kernel (4) may be approxi-mated with any requested accuracy by a finite number ofequations.[36] The definition of a limit cycle in the system (3) and (4)with homogeneous kernels K i,j and F i,j has some sub-tleties discussed in the SM. Supplemental Material: Oscillations in aggregation-shattering processesNumerical versus analytical solutions for the constant kernels
Figure 7 illustrates the application of the fast-integration method for the case of constant aggregation and shatteringkernels ( a = 0). The numerical solution approaches to the steady-state solution which is known analytically [8]: n k = N √ π (1 + λ ) (cid:20) n (1 + λ ) N (cid:21) k Γ( k − )Γ( k + 1) . Here n = λ/ (1 + λ ) is the stationary density of monomers and N = 2 λ/ (1 + 2 λ ) is the stationary density of clusters.The above distribution refers to the case of unit mass density, M = 1, the general case is obtained by multiplying theabove densities by M . For λ (cid:28) k (cid:29)
1, the above exact distribution simplifies to n k = λ √ π e − λ k k − / . (10)Figure 7 demonstrates the high accuracy of the numerical method and the intuitively obvious feature that the smallerthe parameter λ the longer it takes for the system to reach the steady state (when λ = 0, the steady state is neverreached). -12 -8 -4 n k k λ = 0.05 AnalyticalNum. t=400 -12 -8 -4 n k k λ = 0.02 AnalyticalNum. t=2000
FIG. 7: Comparison of the numerical solution with the analytical steady-state solution [8] for the constant kernel K i,j = 1 andmonodisperse initial conditions, n k (0) = Mδ ,k with M = 1, for λ = 0 .
05 (left panel) and λ = 0 .
02 (right panel).
Limit Cycles of kinetic equations with homogeneous kernels
One must be careful while talking about limit cycles for kinetic equations with homogeneous kernels. By definition,a limit cycle is an isolated closed trajectory; this means that its neighboring trajectories are not closed — they spiraleither towards or away from the limit cycle [13]. After this definition, we are usually told that limit cycles can onlyoccur in nonlinear systems; in a linear system exhibiting oscillations closed trajectories are neighbored by other closedtrajectories. We also learn that a stable limit cycle is one which attracts all neighboring trajectories. A system witha stable limit cycle can exhibit self-sustained oscillations [13].Consider a dynamical system that may be written as dn k dt = F k ( n ) , k ≥ n = ( n , n , n , . . . ) and F k ( n ) is given in our case by Eqs. (3) and (4) of the main text. It is importantto note that the reaction terms F k ( n ) are strictly quadratic polynomials for all k , and this fact alone leads to theconclusion that limit cycles in the dynamical system (11) are impossible. Indeed, Eqs. (11) are invariant under thetransformation t → T /M, n k → M N k (12)namely after this transformation Eqs. (11) become dN k dT = F k ( N ) (13)with the same functions F k . Therefore if the dynamical system (11) possesses a limit cycle, we can slightly perturb itby choosing M = 1 + (cid:15) with | (cid:15) | (cid:28) F k ( n ) for all k are homogeneouspolynomials of any degree does not have limit cycles.)We now recall that our dynamical system actually admits an integral of motion, namely the mass density isconserved: (cid:88) j ≥ jn j = 1 . (14)In Eq. (14) we have set the mass density to unity; if the mass density is equal to M , we can make the transformation(12) and then the mass density will be equal to unity.The original dynamical system (11) was considered in the (infinite-dimensional) quadrant R ∞ + = { ( n , n , n , . . . ) | n ≥ , n ≥ , n ≥ , . . . } . (15)But it is more appropriate to reduce (11) to the phase space which is an intersection of (14) and (15). Plugging n = 1 − (cid:80) j ≥ jn j into (11) with k ≥ dn k dt = G k ( n ) , k ≥ . (16)The functions G k ( n ) are still quadratic polynomials, but not strictly quadratic. For example, the term n n turnsinto n − (cid:80) j ≥ jn j n .The dynamical system (16) is defined on ( n , n , . . . ) | n ≥ , n ≥ , . . . ; (cid:88) j ≥ jn j ≤ (17)which is the (infinite-dimensional) simplex. This dynamical system may admit genuine limit cycles. Hence alldynamical systems (11) and (15) with different masses M are equivalent to the generic one, given by Eqs. (16) and(17). In other words, systems with different masses M may be mapped on each other by simple re-scaling. Most ofour simulations have been done for the stepwise initial distribution of the cluster sizes, n j (0) = (cid:40) . ≤ j ≤ j > , (18)with the total mass M = 5 .
5. To illustrate that our limit cycle is unique (up to the numerical precision) we alsoconsider the total mass of M = 3 and mono-disperse initial conditions n k (0) = M δ ,k . In Fig. 6 of the main text it isdemonstrated that the closed trajectories of the n ( t ) − n ( t ) plane coincide for different masses and initial conditionsafter the appropriate re-scaling. This proves numerically the existence of true limit cycle in the system of interest. Analytical approach for the stationary distribution
To find the steady-state distribution of the aggregates sizes, one needs to put ˙ n = ˙ n k = 0 into the left-hand sideof Eqs. (3) and (4) of the main text and solve the following infinite system of algebraic equations: n ∞ (cid:88) i =1 K i n i − λ ∞ (cid:88) i =2 ∞ (cid:88) j =2 ( i + j ) K i,j n i n j − λn ∞ (cid:88) j =2 jK ,j n j = 0 (19)12 k − (cid:88) i =1 K i,k − i n i n k − i − (1 + λ ) n k ∞ (cid:88) i =1 K k,i n i = 0 , k ≥ . We will apply the method of generating functions that has proved its efficiency for similar problems [2, 3, 16]. Namely,we introduce the generating functions C ± a ( z ) and moments M ± a : C ± a ( z ) = ∞ (cid:88) k =1 k ± a n k z k M ± a = ∞ (cid:88) k =1 k ± a n k , (20)Multiplying (19) by z k and summing over all k ≥ C a ( z ) C − a ( z ) + (1 + λ ) zn ( M a + M − a ) − (1 + λ ) ( M a C − a ( z ) + M − a C a ( z )) = 0 . (21)Specializing (21) to z = 1 and taking into account that C ± a (1) = M ± a we obtain M a M − a = 1 + λ λ n ( M a + M − a ) . (22)The tail of the size distribution can be extracted from the asymptotic behavior of the generation functions C a ( z ). Weconsider separately the cases of a < / a > / Kernels with a < / . The tail (10) arising in the context of the model with constant kernel, a = 0, in the casewhen additionally λ (cid:28)
1, suggests that generally steady-state distribution may have a similar tail, n k (cid:39) Ck − τ e − ωk for k (cid:29) , (23)for kernels with a >
0. The amplitudes C and ω and the exponent τ are yet unknown functions of λ and a . Thegeneration functions may be expanded near z (cid:48) → −
0, where z (cid:48) = z/z and z = e ω . One seeks the expansions inthe form [2, 16] C a ( z ) = C a ( z ) + C Γ(1 + a − τ )(1 − z (cid:48) ) τ − a − (24) C − a ( z ) = C − a ( z ) + C Γ(1 − a − τ )(1 − z (cid:48) ) τ + a − , (25)where Γ( x ) is the gamma function and we assume that C a ( z ) = (cid:80) k ≥ k a n k z k < ∞ exists for given a and τ . If wesubstitute the above C ± a ( z ) into Eq. (6) we obtain terms with different powers of the factor (1 − z (cid:48) ). To satisfy thisequation we equate to zero all these terms separately. Then the zero-order terms of (1 − z (cid:48) ) yield C a ( z ) C − a ( z ) − (1 + λ ) ( M a C − a ( z ) + M − a C a ( z )) + (1 + λ ) z n ( M a + M − a ) = 0 . (26)Similarly, the terms of the order (1 − z (cid:48) ) τ ± a − imply the relations C ∓ a ( z )Γ(1 ± a − τ ) − (1 + λ ) M ∓ a Γ(1 ± a − τ ) = 0 . (27)Finally, the rest of the terms should satisfy C Γ(1 + a − τ )Γ(1 − a − τ )(1 − z (cid:48) ) τ − − (1 + λ ) z n ( M a + M − a )(1 − z (cid:48) ) = 0 . (28)This equation is consistent when 2 τ − τ = 32 . (29)The exponent τ is therefore universal (independent on a and λ ). Now we substitute the relations C ∓ a ( z ) = (1 + λ ) M ∓ a , (30)which follow from Eq. (27) into Eq. (26) to yield M a M − a = z λ n ( M a + M − a ) . (31)From Eqs. (31) and (7) then follows, z = e ω = (1 + λ ) (1 + 2 λ ) . (32)The ansatz (23) is expected to work only when λ (cid:28)
1. In this limit (32) gives ω (cid:39) λ − λ + . . . (cid:39) λ , so theamplitude ω is independent on a . Thus the tail of the cluster size distribution reads n k (cid:39) Ck / e − λ k for k (cid:29) . (33)An order-of-magnitude estimate for the constant C may be done as follows. We assume that the distribution (33),which holds true for k (cid:29) k ∼
1, so that ∞ (cid:88) k =1 kn k (cid:39) (cid:90) ∞ Ck / e − λ k dk (cid:39) C √ πλ = 1 , (34)that is, C (cid:39) λ/ √ π ∼ λ . Kernels with a > / . Applying the same analysis as above for a ≥ /
2, one arrives at Eqs. (26)–(28), whichhowever do not lead to consistent results. Indeed, from Eq. (28) it follows that τ = 3 /
2, but C a ( z ) does not exist for a ≥ /
2, so that Eq. (27) may not be satisfied to cancel the terms corresponding to the factor (1 − z (cid:48) ) τ + a − .Although the above approach fails to make consistent asymptotic estimates for a ≥ /
2, our results for a < / a ≥ /
2, the distributionof cluster size has the following form for k (cid:29) n k (cid:39) Ck − / e − λ β k ; it will be used below for the further analysis. Truncating an infinite system of equations by a finite number of equations
The standard problem of numerical solution of Smoluchowski equations is how to approximate an infinite systemof equations by a finite one. When fragmentation is lacking, as in common Smoluchowski equations, the average sizeof aggregates infinitely grows which imposes a principle time limit for the modeled processes. Contrary, in the case ofinterest, the fragmentation of aggregates precludes the formation of very large clusters even for infinitely long time.Therefore the number of equations may be finite. Moreover, using the results for the steady-state distribution, onecan estimate the number of equations needed to describe the system with a given degree of accuracy. Below we show,how the solutions of a formally infinite system (3) of the main text may be adequately represented by these of a finitesystem.Using Eqs. (3) of the main text, we write for the concentration n k ( t ) for 2 ≤ k ≤ L : dn k dt = 12 k − (cid:88) i =1 K i,k − i n i n k − i − (1 + λ ) n k L (cid:88) i =1 K k,i n i − (cid:40) (1 + λ ) n k ∞ (cid:88) i = L +1 K k,i n i (cid:41) . (35)Taking into account that for a ≤ k < L < iK i,k = (cid:18) ik (cid:19) a + (cid:18) ki (cid:19) a ≤ i a + k a < i + i = 2 i and applying the steady-state distribution, n k (cid:39) Ck − e − λ β k ∞ (cid:88) i = L +1 K i,k n i < ∞ (cid:88) i = L +1 in i ∼ (cid:90) ∞ L xCx − e − λ β x dx ∼ CL Erfc( √ λ β L ) √ λ β L .
If the quantity λ β L is large, one can make a further simplification, using the asymptotic relation Erfc( x ) ∼ e − x /x ,which yields, ∞ (cid:88) i = L +1 K i,k n i < C e − λ β L λ β L < ε. (36)Hence, if we choose the number of equations L = N eq ( λ ) such that the above expression is smaller than ε (cid:28) dn k dt = 12 k − (cid:88) i =1 K i,k − i n i n k − i − (1 + λ ) n k N eq (cid:88) i =1 K ki n i , that is, the solution of an infinite system may be approximated with any desired accuracy by the solution of a finitesystem with the appropriately chosen number of equations. In practice, we started with the number of equationsestimated from Eq. (36) for β = 2 and C ∼ λ , as for the steady-state distribution for a < / N eq = 150 ,,