Stochastic gel-shatter cycles in coalescence-fragmentation models
Brennen T. Fagan, Niall J. MacKay, Dmitri O. Pushkin, A. Jamie Wood
eepl draft
Stochastic gel-shatter cycles in coalescence-fragmentation models
Brennen T. Fagan , Niall J. MacKay , Dmitri O. Pushkin and A. Jamie Wood , Department of Mathematics, University of York, York, YO10 5DD, UK Department of Biology, University of York, York, YO10 5DD, UK
PACS – Stability and fragmentation of clusters
PACS – Oscillations, chaos, and bifurcations
PACS – Diffusion and aggregation
Abstract – We describe a new phenomenon in models of coalescence and fragmentation, that of gel-shatter cycles . These are dynamical, unforced, stochastic cycles in which slow, approximatelydeterministic coalescence up to and beyond gelation is followed by abrupt random shattering.We describe their appearance in simulations of stochastic models with multiplicative kernels forcoalescence and spontaneous fragmentation into monomers (‘shattering’). The regime in whichsuch cycles occur is characterized by a cyclicity order parameter, and we provide a simple scalingargument which describes both this regime and those which border it.
It is natural that coalescence of smaller units shouldresult in larger accumulations that may eventually be-come unstable and fragment. A wide range of models forsuch coalescence and fragmentation processes have beencreated, often with overlapping or even contradictory re-sults, in fields as diverse as physical chemistry [1], prob-ability [2], fluid dynamics [3], social grouping [4], bankmergers [5] and terrorist grouping networks [6–8]. Partof this popularity is undoubtedly due to the emergenceof power-law behaviour from simple coalescence, e.g. twoasteroids collide [9], one fish eats another [10] or two in-vestors trade information [11], all of which yield power-law-distributed steady states, i.e. p ( x ) ∝ x − α for α > Fig. 1: Cluster size distribution. A typical time-averaged dis-tribution (black dots) emerging as a result of multiplicative co-alescence and multiplicative shattering, with fragmentation re-stricted to clusters larger than 10 . The total mass is M = 10 .(a) Complementary cumulative distribution of clusters, i.e. thenumber of clusters larger than a certain size. (b) Mean clus-ter density. This follows a truncated α = power-law (bluedashed line). We also show the distribution at a single timewhose power-law exponent is closest to the median (red line). p-1 a r X i v : . [ m a t h . D S ] F e b agan et al. dn k dt = 12 k − (cid:88) i =1 K ( i, ( k − i )) n i n k − i − n k ∞ (cid:88) i =1 K ( i, k ) n i − F ( k ) n k + ∞ (cid:88) i = k +1 F ( i ) b ( i, k ) ik n i . (1)in the absence of cluster influx/outflux [17]. Most recently,it has been shown that deterministic temporal oscillationscan arise in a class of coalescence-fragmentation (C-F) pro-cesses where clusters grow or shrink by addition or deletionof monomers [18]. All of these findings predicted determin-istic oscillations and relied on mean-field descriptions.Here we use simulations to show that stochastic cyclicaldynamics – in the form of steady accumulation followedby sudden shattering – emerge in an important class of C-F processes. We consider spontaneous fragmentation, notinduced by collision. Such C-F processes have been usedto study polymerisation/depolymerisation reactions [19],cloud droplets spectra [20], and assembly/disassembly ofcell membrane microdomains [21], among other applica-tions. We focus on a particular example of size-biasedcoalescence and fragmentation processes which is widelyused in applications. We provide simple scaling argumentsto explain the observed dependence of the recurrence timeon the system size.Mathematical descriptions of the problem can begrouped into two broad classes: mean field approachesand stochastic models (including stochastic simulations).Coalescence processes, also known as coagulation or aggre-gation processes, were introduced to mathematics a cen-tury ago [22], while the stochastic equivalent, the Marcus-Lushnikov process, arose in the 1960s [23, 24].First, for reference and convenience, we present an adap-tation of Smoluchowski’s equations for coalescence andfragmentation, modified from reference [25]. In equation(1), t is time and n k is the mean field density of clustersof (discrete) size k . Larger clusters are generated whentwo smaller clusters of sizes i and j randomly coalesce, inproportion to a kernel K ( i, j ), to form a cluster of size i + j . Clusters are reduced in size by a separate processof fragmentation, which occurs with the rate F ( i ) for acluster of size i . The function b ( i, j ) describes what pro-portion of mass of the cluster of size i goes to a clusterof size j ; it is normalised so that (cid:80) i − j =1 b ( i, j ) = 1. Thiscondition guarantees that the mass is conserved in a frag-mentation event. Since here fragmentation is spontaneousrather than collision-induced, the terms describing frag-mentation are linear in the density of clusters.The mean field approaches provide the basic languagefor cluster-size distribution modelling but ignore finite sys-tem size. However, in any computational implementationfinite-size effects associated with the total system size M will be important. In particular, the dynamics of pure co-alescence, when the coalescence rate grows sufficiently fastas a function of the sizes of the coalescing clusters, leads toa giant cluster of a size comparable to M , a gel [2, 26–28].[27]. In the thermodynamic limit, when M → ∞ and the gel size is infinite, gelation goes beyond mean field theory(1) and stochastic modelling is essential for capturing theensuing dynamics.A particularly important coalescence kernel allowing gelformation naturally emerges in describing random networkgrowth. When two nodes are randomly connected by anedge per unit time, two clusters of connected nodes coa-lesce at the rate K = ˆ K ( i/M )( j/M ), where ˆ K is constant,the system size M is the total number of nodes, and i/M and j/M are the probabilities of picking clusters of sizes i and j from the population. This kernel is often referredto as the multiplicative coalescence kernel.Similarly, the multiplicative fragmentation rate F ( i ) isproportional to the cluster size i . The particular case inwhich the fragmented clusters disaggregate into monomersis referred to as shattering [13, 14] (We adopt the termi-nology used by reference [13], but note that shattering hasbeen used to describe alternative phenomena, e.g. [29].)This form of fragmentation can be conveniently simulatednumerically by picking a random node from the popu-lation of all nodes and shattering the cluster it belongsto, removing all edges amongst the formerly connected Fig. 2: Stochastic coalescence-fragmentation cycles. (a) Num-ber of clusters N versus maximum cluster size k max . Thecounterclockwise trajectory shows periods of relatively gradualgrowth of the largest cluster followed by its abrupt shatter-ing. (b) The time-dependent power-law exponent α obtainedby fitting the instantaneous cluster-size distribution to a (trun-cated) power-law using maximum likelihood estimation. Thecoloured region corresponds to the trajectories shown in (a).(c) Maximum cluster size as a function of time. Here ˆ K = 0 . F = 0 .
01 and the total populations is M = 10 . An extendedand animated version of this figure is included with the Sup-plemental Material. p-2el-shatter cyclesnodes. The corresponding fragmentation rule is given by F ( i ) = ˆ F ( i/M ) , b ( i, j ) = δ j , where ˆ F is constant and i/M is the probability of picking a cluster of size i .Simple kernels, especially those that are multiplicative,often result in solvable steady state systems, and the aboveforms are no different. Under the assumption of a steadystate, dn k /dt = 0 for all k , set ρ k = kn k /M so that for k ≥ ρ k = ˆ K ˆ F + ˆ K (cid:80) ∞ i =1 ρ i k − (cid:88) i =1 ρ i ρ k − i . Provided that (cid:80) ∞ i =1 ρ i is well-behaved, which empiricallyis the case when ˆ F is sufficiently large compared to ˆ K ,we have a variation of the standard recurrence relation ofthe Catalan numbers [30, 31]. One standard method ofsolving such a system is to use generating functions andseries expansion [7, 11]. Writing the result in terms of theCatalan numbers C n and writing the prefactor as γ , wehave ρ k = C k − γ k − ρ k . (2)Stirling’s law then yields the standard result of a trun-cated power-law with exponent − /
2. The ρ term canbe evaluated by considering the total amount of mass inthe system to be fixed. Many other, related solutions existfor the steady-state distribution, dependent on the preciseformulation of Equation 1 [3, 12, 13, 15].The typical stochastic dynamics and its emerging cyclesare seen in Figure 2(a), which shows a sample of the trajec-tory in the ( k max , N ) plane, where N is the total numberof clusters and k max is the size of the largest cluster. Eachdistinct counterclockwise cycle begins in the top-left witha predominance of monomers. A period of gradual growthof k max is accompanied by a decrease in N , taking the tra-jectory slowly down and to the right before acceleratingrightward through gelation until growth ends in an abruptrandom jump back to top-left with ∆ N ≈ − ∆ k max . Thisclearly evidences that 1) coalescence leads to the forma-tion of very large clusters during the periods of gradualcluster size growth and 2) such periods end in shatter-ing of the largest cluster. However, the cycling dynam-ics are not limited to the growth and shattering of thelargest cluster but involves the whole cluster size distri-bution. Indeed, our simulations show that the cluster sizedistribution is very broad and can be fitted to (truncated)power-laws with the time-dependent exponent α ( t ) cyclingin the range 2 . − .
9. During the period of coalescence,the cluster size distribution is broadening and α decreasesappreciably.To the best of our knowledge, such stochastic gel-shattercycles have not been previously described. We hypothe-size that they should be a generic feature of C-F systemswith gel-forming coalescence and sufficiently strong frag-mentation.The gel-shatter cycles are strongly dependent on thenumber of monomers at complete disaggregation, whichwe term the system size M . Figure 3 presents a ‘heat map’ Fig. 3: Cluster summary statistics. System sizes are M =3 × and 3 × (left to right), and fragmentation ratesˆ F = 10 − and 10 − (top to bottom, with coalescence ratesˆ K = 1 − ˆ F ). The number of clusters N and the maximum clus-ter size k max are normalised by M . Darker regions representstates that emerge more often. The red border line denotes theboundary of the region visited by the system. A more detailedfigure is present in the Supplemental Material. of the times the simulation spends in different regions ofthe ( k max , N ) plane (normalised by M for comparison be-tween different populations). For low fragmentation ratesand small system sizes (top-left corner) the system visits abroad shoulder-like region. The ‘elbow’, the change in gra-dient of its lower boundary, corresponds roughly to gela-tion. As fragmentation rate and system size increase, thebroad distribution collapses to a very small region withinwhich the stochasticity is not visible as distinct cycles.The mean recurrence time (cid:104) t r (cid:105) is defined to be the num-ber of computational steps between successive shatteringsof the largest cluster, averaged across simulations with thesame parameters. This can be thought of as the averageduration of a cycle in the gel-shatter regime. The depen-dence of the mean recurrence time (cid:104) t r (cid:105) on system size M for varying fragmentation rates ˆ F shows crossovers be-tween three distinct regimes: for small systems and smallfragmentation rates, (cid:104) t r (cid:105) ∼ M while modestly increas-ing M or ˆ F yields (cid:104) t r (cid:105) ∼ M / , before crossing over into aregime with superlinear growth of (cid:104) t r (cid:105) with M . It is a pri-ori possible that these regimes are due to finite size effects.In our simulations, finite-size effects are strong and appeareven for systems of size M = 10 . Furthermore, there aremultiple time scales present, and the importance of thesecould vary across rates and system sizes. These curvesshow data collapse when ˆ F (cid:104) t r (cid:105) is plotted as a function ofthe dimensionless parameter r = ˆ F M ˆ K . (3)p-3agan et al.
Fig. 4: Observed data collapse in terms of r . Points are plotted for M = 10 , , and 10 , but data were additionallygathered at M = 3 × , × and 3 × . (a) When plotting the mean recurrence time (cid:104) t r (cid:105) multiplied by the fragmentationrate ˆ F , we see that the middle region of unforced gel-shatter cycles has strong data collapse following a nearly linear trend onlog-log axes. (b) When plotting the order parameter K , we observe a plateau of K in the region of unforced gel-shatter cycles. The physical meaning of r is the ratio of the characteristictimes of gelation T g ∼ M/ ˆ K and shattering T f ∼ / ˆ F .Figure 4(a) demonstrates that (cid:104) t r (cid:105) = ˆ F − g ( r ) , (4)where the scaling function g ( r ) ∼ r (cid:28) .
1. For largervalues of r , the scaling function crosses over to g ( r ) ∼ r / and this behavior persists for almost four orders ofmagnitude of r . Finally, for r (cid:38) , the scaling breaksdown.In order to explore further the nature of the stochasticgel-shatter cycles and distinguish cyclical dynamics fromacyclical stochastic fluctuations, we introduce a cyclicityorder parameter, K , defined to be the number of computa-tional steps that result in growth of the maximum clustersize k max minus those that result in reduction, divided bythe total number of computational steps. Hence, K mustlie in the range − K , in contrast, indicatesmany steps of growth followed by abrupt collapse, therebycharacterizing gel-shatter cycles. (Large negative K wouldindicate many fragmentation events followed by rare co-alescence. This is possible only if fragmentation eventsare frequent and small, and we do not expect to observeit with our chosen ˆ K and shattering fragmentation ker-nel F .) The order parameter K depends non-trivially on r , see Figure 4(b). Here K ( r ) forms hump-shaped curveswith the maximum around r = 10. The order parameteris appreciably large, K > .
1, for 0 . < r < , whichconstitutes the middle regime with (cid:104) t r (cid:105) ∼ M / .We can now identify the different scalings as three phys-ical regimes: (i) Weak fragmentation or forced cycles.
Here T g (cid:28) T f and coalescence quickly results in a singlegel cluster. The gel is shattered at the rate ˆ F andhence the emerging cycles have the mean recurrencetime (cid:104) t r (cid:105) ∼ ˆ F − . This regime is physically unsur-prising.(ii) Unforced gel-shatter cycles.
Here T g ∼ T f andthe cycles arise from a non-trivial interplay betweenthe dynamics of gelation and shattering. The cyclic-ity K reaches a maximum in this regime. The pres-ence and extent of this regime in gelling C-F systemsis our main finding.(iii) Fragmentation-dominance.
Here T g (cid:29) T f andfragmentation is so strong as to preclude formationof large clusters. This annihilates the gel-shatteringcycles, leaving only stochastic variation in a smallregion.We now provide a simple scaling argument to explainthe behaviour of the mean recurrence times and gain fur-ther insight into the gel-shatter cycles.Assume the dynamics are dominated by continuousgrowth and shattering of the largest cluster. Writing thesize of the largest cluster k max at time t as m ( t ), we trackthe largest cluster and the probability, dependent on m ( t ),that it shatters at a given time. The probability P ( t ) thatthe largest cluster is not shattered by time t is governedby dP ( t ) dt = − ˆ F M − m ( t ) P. (5)Hence the probability density for shattering at time t is p ( t ) = − dP ( t ) dt = ˆ F M − m ( t ) e − ˆ F M − (cid:82) t m ( τ ) dτ . (6)p-4el-shatter cyclesIn the forced-cycles regime (i), coalescence quickly leadsto a single large cluster of size M , hence m ( t ) = M and p ( t ) = ˆ F exp (cid:16) − ˆ F t (cid:17) . This is the exponential distribution with mean (cid:104) t r (cid:105) = ˆ F − ;thus g ( r ) ∼ , r (cid:28) m ( t ) consistent with dimensional argumentsis linearity. We assume that m ( t ) = c ˆ Kt , where c > p ( t ) = c ˆ F ˆ KM − te − c ˆ F ˆ KM − t / . This is the Rayleigh distribution with scale parameter (cid:16)(cid:112) c ˆ F ˆ KM − (cid:17) − and mean (cid:104) t r (cid:105) = (cid:113) πM/ (2 c ˆ F ˆ K ) . (7)Hence, g ( r ) ∼ r / in regime (ii). Also, it is straightfor-ward to show that the largest realisable cluster size scalesas M r − / . This scaling explains the shortening of the‘elbow’ on Figure 3 for larger M and ˆ F .This work reports a distinctive new phenomenon, thatof stochastic gel-shatter cycles . Gel-shatter cycles areobserved explicitly in simulations of multiplicative (size-biased) coalescence and spontaneous-shattering fragmen-tation kernels when the time scales for these are in balance.We expect gel-shatter cycles to be a general phenomenonbeyond that of multiplicative kernels. The prerequisitequalities appear to be three-fold: gelation, which promotesthe growth of a single cluster that dominates the system; astrong form of fragmentation, which can reset the systemto a pre-gel state; and balanced time scales, so that the gelis neither instantly removed nor aggregates the entire sys-tem. Gelation is a broad and well-studied phenomenon,occurring with many coalescence kernels. Shattering isperhaps the strongest assumption we make, but we antic-ipate that this can be relaxed to a sufficiently strong formof mixed fragmentation. Finally, in any system studied forits tension between coalescence and fragmentation (as op-posed to just one of these alone), the time scales of theseare likely to be comparable. In such circumstances prac-titioners should be wary of presuming the existence of asteady state, and conscious that gel-shatter cycles may bepresent.On this basis we anticipate that gel-shatter cyclesare a ubiquitous phenomenon. This salient novel fea-ture emerges from a range of simulation studies origi-nally intended to explore the robustness of the Smolu-chowski coalescence-fragmentation model and its variants(see Supplemental Material). Whilst our conclusion is thatmodels of this type are robust to perturbations of vari-ous kinds – confirming the common view [6–8, 12] – wenote the ubiquitous presence of unforced stochastic gel-shatter cycles in finite-size systems. This is relevant tothe application of these models to physical and social phe-nomena, where the gel-shatter cyclicity may be a natural dynamical feature, as well as affecting the robustness ofthe computationally-fitted exponents. The presence of anon-trivial dynamic underlying a presumed steady stateneeds to be considered when models of this type are usedin applications, and is an intriguing avenue for future the-oretical research. ∗ ∗ ∗ We would like to thank Stephen Connor and GustavDelius for interesting discussions. We are also grateful totwo anonymous referees for their constructive comments.
REFERENCES[1]
Spicer P. T. and
Pratsinis S. E. , AIChE Journal , (1996) 1612. https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.690420612 [2] Aldous D. J. , Bernoulli , (1999) 3. [3] Pushkin D. O. and
Aref H. , Physics of Fluids , (2002) 694. https://doi.org/10.1063/1.1430440 [4] Johnson N., Spagat M., Restrepo J., BohorquezJ., Suarez N., Restrepo E. and
Zarama R. , ArXivPhysics e-prints , (2005) .[5]
Pushkin D. O. and
Aref H. , Physica A: StatisticalMechanics and its Applications , (2004) 571 . [6] Bohorquez J. C., Gourley S., Dixon A. R., SpagatM. and
Johnson N. F. , Nature , (2009) 911. [7]
Clauset A. and
Wiegel F. W. , Journal of Conflict Res-olution , (2010) 179 . http://dx.doi.org/10.1177/0022002709352452 [8] Johnson N. F., Zheng M., Vorobyeva Y., GabrielA., Qi H., Velasquez N., Manrique P., Johnson D.,Restrepo E., Song C. and
Wuchty S. , Science , (2016) 1459. http://science.sciencemag.org/content/352/6292/1459 [9] Tanaka H., Inaba S. and
Nakazawa K. , Icarus , (1996) 450 . [10] Datta S., Delius G. W., Law R. and
Plank M. J. , Journal of Mathematical Biology , (2011) 779. https://doi.org/10.1007/s00285-010-0387-z [11] D’Hulst R. and
Rodgers G. J. , International journalof theoretical and applied finance , (2000) 609 .[12] Ruszczycki B., Burnett B., Zhao Z. and
JohnsonN. F. , The European Physical Journal B , (2009) 289. http://dx.doi.org/10.1140/epjb/e2009-00354-5 [13] Birnstiel T., Ormel C. W. and
Dullemond C. P. , Astronomy & Astrophysics , (2011) A11. https://doi.org/10.1051/0004-6361/201015228 [14] Connaughton C., Dutta A., Rajesh R., SiddharthN. and
Zaboronski O. , Phys. Rev. E , (2018) 022137. p-5agan et al. https://link.aps.org/doi/10.1103/PhysRevE.97.022137 [15] Kyprianou A. E., Pagett S. W. and
Rogers T. , An-nales de l’Institut Henri Poincar´e, Probabilit´es et Statis-tiques , (2018) 1134.[16] Ball R. C., Connaughton C., Jones P. P., RajeshR. and
Zaboronski O. , Physical Review Letters , (2012) . https://doi.org/10.1103%2Fphysrevlett.109.168304 [17] Matveev S. A., Krapivsky P. L., Smirnov A. P.,Tyrtyshnikov E. E. and
Brilliantov N. V. , Phys.Rev. Lett. , (2017) 260601. https://link.aps.org/doi/10.1103/PhysRevLett.119.260601 [18] Pego R. L. and
Vel´azquez J. J. L. , Nonlinearity , (2020) 1812. https://doi.org/10.1088%2F1361-6544%2Fab6815 [19] Blatz P. and
Tobolsky A. , The Journal of PhysicalChemistry , (1945) 77.[20] Pruppacher H. R. and
Klett J. D. , Microphysics ofClouds and Precipitation (Springer) 2010.[21]
Richardson G., Cummings L. J., Harris H. and
O’Shea P. , Biophysical Journal , (2007) 4145.[22] von Smoluchowski M. , Zeitschrift fuer physikalischeChemie , (1916) 129 . http://publikationen.ub.uni-frankfurt.de/frontdoor/index/index/docId/13699 [23] Marcus A. H. , Technometrics , (1968) 133. [24] Lushnikov A. A. , Journal of Colloid and InterfaceScience , (1978) 276 . [25] Spouge J. L. , Mathematical Proceedings of the Cam-bridge Philosophical Society , (1984) 351–357.[26] Wattis J. A. D. , Physica D: Nonlinear Phenomena , (2006) 1.[27] Ziff R. M., Hendriks E. M. and
Ernst M. H. , Phys.Rev. Lett. , (1982) 593. https://link.aps.org/doi/10.1103/PhysRevLett.49.593 [28] Lushnikov A. , Physica D: Nonlinear Phenomena , (2006) 37 .[29] McGrady E. D. and
Ziff R. M. , Phys. Rev. Lett. , (1987) 892. https://link.aps.org/doi/10.1103/PhysRevLett.58.892 [30] Hilton P. and
Pedersen J. , The Mathematical Intelli-gencer , (1991) 64 . https://doi.org/10.1007/BF03024089 [31] OEIS Foundation Inc , Catalan numbers (2019). https://oeis.org/A000108https://oeis.org/A000108