Universality of Cluster Dynamics
Abstract
We have studied the kinetics of cluster formation for dynamical systems of dimensions up to interacting through elastic collisions or coalescence. These systems could serve as possible models for gas kinetics, polymerization and selfassembly. In the case of elastic collisions, we found that the cluster size probability distribution undergoes a phase transition at a critical time which can be predicted from the average time between collisions. This enables forecasting of rare events based on limited statistical sampling of the collision dynamics over short time windows. The analysis was extended to Lnormed spaces () to allow for some amount of interpenetration or volume exclusion. The results for the elastic collisions are consistent with previously published lowdimensional results in that a power law is observed for the empirical cluster size distribution at the critical time. We found that the same power law also exists for all dimensions , 2D L norms, and even for coalescing collisions in 2D. This broad universality in behavior may be indicative of a more fundamental process governing the growth of clusters.
pacs:
05.20.Dd, 05.65.+b, 45.70.Vn, 45.50.Tn, 45.70.Vn, 89.75.kI Introduction
This paper is a study of the statistical behavior of the dynamics of clusters which are allowed to interact through elastic collisions or by coalescence. The elastic collision dynamics are based on a ballistic billiard model analyzed theoretically by Sinai bib:sinai (); bib:sinai2 (). Cluster growth and selfassembly processes are relevant to a variety of research fields, including chemistry, materials science, physics and earth sciences. The study of such random processes can reveal information on the nature of collective interactions as well as make predictions on the occurrence of rare and catastrophic events. Early theoretical studies on cluster dynamics originate in the work of Bogoliubov, who showed that in the gas phase, groups of particles with shortranged interactions behave like independent clusters bib:bogolyubov (). Sinai provided a proof of cluster dynamics for colliding billiards for one dimensional (1D) systems bib:sinai (), and subsequently for higher dimensions (restricted to sufficiently low densities). Sinai also proved ergodicity of the classical billiard model bib:sinai2 (). The statistical properties of cluster dynamics has been studied for the 2D case with frictionless elastic billiards gabriel (). In this paper we extended this statistical analysis of ballistic billiards to higher dimensions () up to , higher densities (), Lnormed distance metrics () and to the case of coalescing billiards. We have found a high degree of universality which suggests that the dynamics of clustering are relatively independent of the details.
In an ensemble of interacting particles we may observe a phase transition where a dominant cluster emerges gabriel (). In a classical Sinai billiard consisting of elastic collisions, the phase transition in the empirical density of clusters is not necessarily associated with a phase transition of the physical system in the traditional sense, such as a transition from liquid to solid or gas to liquid as function of temperature or pressure. Instead, one observes a change in the empirical cluster density – which plays the role of the order parameter – as function of time. Thus, it is indicative of the dynamics of the motion rather than a configurational change resulting from the variation of an intensive variable. The collisions between billiards represent the interactions between parts of a system, and the transition to a dominant cluster that emerges is a manifestation of the events leading to a major catastrophic event. In the case of coalescence, the phase transition in the probability density can be associated with a physical change in the properties of the system. Recent examples of the analysis of phase transitions in probability densities include earthquake prediction keilis (); gabriel2 (); gabriel3 (); zal3 (); zal4 (), economic modeling keilis2 () and models of river networks zaliapin (). Established premonitory patterns have allowed the modeling of events in complex systems to be predicted using observed background activity gab2 (). The prospect of predicting or controlling critical events occurring in a dynamic and complex environment is of broad interest.
In the first part of the paper, we expand the study of dynamical phase transitions to higher dimensional elastic billiards and describe the statistics of the collisions in dimensional Lnormed spaces, with and . The first result which emerges is the existence of the phase transition in dimensions greater than =2, for higher densities, and for different Lnormed spaces (). In the Euclidean norm case of Sinai billiards, the critical time appears to be independent of the dimensionality of the system (). Instead, this critical time solely depends on the average time between collisions . Another notable finding is that the empirical cluster distribution at the critical point obeys a power law across all dimensions, densities and norms with the same exponent.
In the second part of the paper, we allow the billiards to coalesce and form larger clusters. These collisions can be analyzed using a binary tree model zaliapin () first developed to analyze environmental transport in river networks. Coalescence and coagulation are phenomenon that are present in many areas of chemistry kang (); ghosh (). Theories of coalescence date back to the work of Smoluchowski smoluchowski (); smoluchowski2 () in the early 20th century, establishing the evolution of the concentration, , of clusters of mass using a master equation of the form (discrete case):
(1) 
where is the interaction kernel, which is dependent on the collision process of mers and mers. The first term predicts an increase in due to coalescence of an mer and mer; the second term deals with the decrease in due to mers coalescing with clusters of different sizes kang (). The theory is based on two important assumptions: coalescence upon collision and the absence of hydrodynamic interaction between the different mers wang (). Recent work in the field has led to corrections to Smoluchowski’s equation. Such works include film drainage theory chesters () and studies considering a hydrodynamic interaction term wang (); zeichner ().
The coalescence process we analyze is similar to that of Smoluchowski in the sense that there are no interparticle interactions except for coalescence events which occur upon contact, and the process begins at with a monodisperse collection of monomers. The model is found to exhibit similar properties with regards to universality of the phase transition as the elastic model. The results suggest that coalescing processes are governed by principles similar to that of noncoalescing billiards.
Ii Model and Definitions
We start with a Sinai billiard gabriel () in dimensions, involving spheres positioned inside a frictionless hypercubic domain. The total mass, , of the billiards in each case is 1.0
with radius . The domain is the set of points
(2) 
The density of billiards within this domain is:
(3) 
where is the volume of a dimensional hypersphere of radius and is the volume of the hypercubic domain.
Clusters are defined in the elastic collision model using the notion of a cluster sinai2 (). Time is measured by the variable , and specifies an interval of time. A cluster is a group of billiards which have effected each others’ kinematics in the previous time interval . A neighbor is defined as two billiards at time which have collided during the time interval . The set of all billiards which have interacted and are linked by neighbor relationships, are called a cluster. A cluster’s mass, , is the sum of the mass of the billiards which make up the cluster. We use the notation at time to represent the mass of the largest cluster (by mass) and as the total number of clusters at the time .
For the coalescing case, clusters are defined as neighbors or clusters only if they coalesce upon collision. Each cluster corresponds to one billiard in the hypercubic domain. is the mass of the largest cluster at time and is the total number of clusters at time .
ii.1 Kinetics of Collisions
We use a ballistic colliding billiard model for both the elastic and coalescing cases. We use an a priori method to detect collisions in the Euclidean norm case and coalescing case. When varying the norm (), we use an a posteriori method to detect collisions. Each of these are further described below.
ii.1.1 Sinai Billiard
Elastic hard sphere collisions are enforced for the Sinai billiard. Total energy, , and momentum, , of the billiards in the system remain constant:
(4) 
where . Incident and reflective angles are identical for billiardwall collisions.
In constructing the simulation, an a priori method was used to calculate the next time, , of collision between two billiards or between billiard and wall. The system was advanced to this time, , the velocities redefined for colliding billiards, clusters recorded, and the next collision time, , computed.
ii.1.2 LNormed Spaces
We also investigated the effects of using an Lnorm (). The motivation for this is to consider particles that are of irregular shape and surfaces not unique to Euclidean space in collisions. For the application of Lnorms in this context, see gamba (). The L norm is
(5) 
When dimensional billiards collide in any Lnorm, momentum and energy are conserved by reassigning velocities as follows:

The vector normal to the sphere’s surface at the point of collision is calculated between two colliding billiards, and . First, by calculating the normal component along each dimension, using the billiard positions and .
(6) where . We then calculate the associated unit normal vector :
(7) 
At the time of collision, the initial (i) velocity of is given by . We denote the velocity for analogously. The initial relative velocity of and is calculated and dotted with the unit vector, , to find the speed, associated with the impulse.
(8) 
The impulse, , is
(9) 
Velocities are reassigned for the two billiards to find the final velocity (f) after collision.
(10)
In implementing the L norm dynamics, an a posteriori method is used where a constant time interval, , was used to advance the billiards at each time step. Billiardbilliard intersections or billiardboundary overlaps are used to identify collisions. When those collisions are identified, velocities are redefined using the method described above. The model then proceeds to the next time step. The model was verified against the a priori model in the (Euclidean) case in order to find a time step that is appropriate for simulation speed.
ii.1.3 Coalescing Billiard
Coalescing billiards were introduced to model polymerization reactions and selfassembly processes. In a different but mathematically similar context, Zaliapin zaliapin () studied transport in river networks. This description of coalescing processes can include, for example, emulsions of oil in water, or colloidal particles flocculation in percolation analysis hasmy (). As a motivation for the analysis one may consider a growing spherical polymer or colloid particle. The probability of coalescence for collisions of particles is assumed to follow an Arrhenius law:
(11) 
where is a temperature and is the activation energy required for coalescence. For individual collisions is the conditional probability of a single coalescence event involving two billiards in a collision taking place where is fixed and is the total kinetic energy of the two colliding billiards. Comparing a uniformly distributed random number in the interval to , coalescence of the two spheres proceeds if ; otherwise, an elastic collision occurs. In the latter case, the kinematics of the elastic collision are the same as that described for the Sinai billiards. In the former case (event of coalescence), the two coalescing billiards, designated daughters and , join to form one billiard, the parent, . In doing so, the radius, , of the parent is defined in terms of the volumes of both initial clusters as to maintain constant density of billiards in the domain of the dimensional system,
(12) 
Mass is also conserved, with , where and are the masses of daughters and respectively. The resulting magnitude of velocity of the parent billiard is defined by conservation of energy:
(13) 
The direction of the parent’s velocity vector is defined as that which results from a completely inelastic collision between the two daughter billiards where momentum is conserved. The parent’s center, , is defined as the center of mass of the two daughter billiards. If the parent extends beyond the boundary of the system upon definition, its center is redefined perpendicular to the boundary edge so that it is completely within the confines of the system. In the event that the parent overlaps with another billiard in the system, that billiard undergoes a 1.0 probability collision with the parent in the same manner as described above.
As time evolves and billiards collide and coalesce, we obtain a binary tree of coalesced billiards. Each cluster is an individual growing billiard on the surface. The total number of clusters is designated as at time and represents the mass of the largest cluster. Due to the computationally demanding nature of this simulation, we investigated only the Euclidean distance metric (L) together with the a priori method.
ii.2 Model Parameters
Both models were simulated varying , the number of billiards ; , the density ; , the activation energy satisfies . At time , nonoverlapping billiards are randomly placed in the volume
(14) 
Particles are assigned an initial Maxwellian distribution:
(15) 
Several L norms in the range were investigated. The temperature is held constant at .
ii.3 Phase Transition in Cluster Dynamics
During an interval the cluster distribution evolves continuously and without gaps. In this regime, the size of the largest cluster at any time, , is not substantially larger than that of the second largest cluster, . A dramatic change (Figure 1) occurs at typical of a phase transition in which a dominant cluster appears. The following definition of the critical time has been proposed gabriel ()
(16) 
This definition has the operational disadvantage of not being a stopping time, meaning that the occurrence of the phase transition cannot be decided based on previous knowledge of history available until the present time. Our results suggest that despite the lack of a stopping time definition, perhaps an even more important observation is the empirical dependence of on the average time between collisions, , highlighting a potential way to predict its occurrence given some amount of statistical sampling. The prediction of is important for the analysis of earthquake events keilis (); gabriel2 (); gabriel3 (); zal3 (); zal4 (), interacting billiards sinai (), economic models keilis2 () and geographical river models zaliapin (). In each of these instances, it has been suggested that cluster dynamics could be used to further understand the system or identify and predict the occurrence of a critical event.
Iii Results
iii.1 Onset of the Phase Transition
In Fig. 1 we vary the dimension with initial parameters , . For each case, D, a similar characteristic formation of clusters occurs with respect to time. The 2D case is in agreement with results previously published by Gabrielov et al. gabriel (). For each a phase transition is found at a critical time, , associated with a rapid increase in the growth of a dominant cluster, .
As the dimension is varied, while keeping other parameters constant, the time scale for the evolution of cluster dynamics is altered, as seen in Fig. 1. There are two contributions leading to this change.

High packing densities become less available for packing spheres at higher dimensions. This point is discussed by Skoge et al. skoge ().

The critical time, , is found to be dependent only on the average time between billiard collisions, for 2. This empirical dependence is shown in Fig. 3 for approximately equal at constant .
The average time between collisions, is calculated from
(17) 
where is the total number of billiardbilliard collisions that have taken place in the time interval .
For fixed density the average time between collisions first increases, then decreases with increasing , as explained by the first point. Due to the critical time’s dependence on , it scales accordingly. There is a general trend that the dimension, , for which reaches its maximum value, increases with decreasing density. Thus, we observe similar cluster distributions at each dimension, scaled by a factor. From our knowledge of the behavior of systems of constant density and varying dimension, we have that as the density decreases, the largest critical time is associated with higher dimensions .
These empirical results provide a predictive tool to estimate the critical time of the system: by estimating from limited statistical sampling, one may predict by using a multiplicative factor. This multiplicative factor is computed from our data in Fig. 4. Data collected and fit found that for = 1000, = 2, the critical time, , as a function of average time between collisions, , for 2D data is = 0.63 + 471.5, for 4D data is = 0.24 + 485.22, and 6D data is = 0.200 + 454.40. This relation for must be calculated for each value of . This gives an expectation value for the critical time, based on parameters of the system. The error for actual measurements from the expectation value increases in magnitude proportional to the critical time. The fractional uncertainty for each data point, :
(18) 
was calculated. The fractional uncertainties of the points form a normal distribution with a standard deviation of for the 2D case, for the 4D case, and for the 6D case. Accounting for this error in the fit relationships, the values fall within each other’s standard deviation for points above . We performed additional runs with for 2D, 4D, and 6D and found that the previous fits and standard deviations of the fractional uncertainties , from the lower data represented well this higher data. In addition, we considered the standard error, comparing the mean critical time for sets of ten runs with equal density to the calculated standard error in time between collisions for that set, . As the critical time increases, so does the standard error of , as we expect from the broadening data as the critical time increases (Fig. 4).
Another method for identifying the phase transition, which was proposed in gabriel (); zaliapin (), analyzes the distribution of clusters to as function of time . In Fig. 2 we plot the empirical cluster density distribution, using fifty (50) independent trials with identical initial parameters and compare the magnitude of each cluster () at a time with the empirical probability density of cluster size. This empirical density for the cluster size distribution evolves over time according to a power law with exponential taper gabriel ():
(19) 
where the critical exponent is , and is a timevarying function expressing the tendency for exponential tail in the sample. As approaches , diverges and a pure power law emerges. The powerlaw fit to our results is shown by a straight line in Figure 2. The premonition of catastrophic events could be done by sampling this statistical distribution over time and ascertain the approach to a power law. This provides an indicator of the imminence of the phase transition.
Nearly identical empirical cluster distributions are found for each dimension at the critical time (Fig. 2). As for the 2D case of gabriel (), varying the density did not result in changes to this distribution. For each dimension, the cluster distribution evolves according to the power law with exponential tail until the critical time, (see Fig. 5). The unique distribution at the critical time indicates a critical exponent for the system at the time of phase transition. The critical exponent is consistently in each dimension studied for elastic collisions.
The usefulness of this premonitory sign in predicting the imminence of a phase transition is determined by the speed at which the power law is approached. We do this by tracking clusters over time and taking the logarithm of both cluster size and density, and approximate the form of the cluster distribution (Eq. 19) as a low order polynomial in at some time :
(20) 
The timeevolution of the cluster distribution is recorded according to Eq. (20), fixing , and fitting the parameters for each . A pure power law emerges at the critical time, as 0.
Cluster density distributions were created by running simulations fifty times for various parameters as reported in Table 1. Results for the timeevolution of the coefficient are presented in Fig. 6. The parameter is seen to decay approximately exponentially. The point, denoted , can be used as a predictive criterion. Fitting the values using , we find values for tabulated in Table 1. The data fits in Table 1 show that a pure power law fit is first approached approximately two to four standard deviations, , prior to the mean, , for fifty runs.
at  

2  0.1  1000  1.271  0.1755  0.1324  0.16601  0.84 
2  0.01  1000  3.1584  0.3167  0.1525  0.410102  2.3 
2  0.001  1000  7.8147  1.0354  0.1312  1.3766  6.0 
2  0.01  2500  3.258  0.2491  0.153  0.351295  2.47 
4  0.1  1000  0.70751  0.09624  0.1339  0.126561  0.54 
4  0.01  1000  5.6819  0.8655  0.1287  0.74666  4.18 
4  0.001  1000  32.208  3.931  0.1146  4.77281  26.6 
4  0.01  2500  7.5752  0.6641  0.1306  0.680786  6.175 
6  0.01  2500  1.9977  0.2280  0.125  0.266022  1.54 
6  0.01  1000  1.3255  0.2133  0.1147  0.203752  1.02 
6  0.001  1000  11.708  1.516  0.1242  1.86351  9.6 
Future research may consider the possible existence for an upper critical dimension where the is no longer valid. Previous work has found for a particular system of coalescing clusters that there is an exact solution of the Smoluchowski equation valid for kang (). It is possible that the reason for the observed universality for could be that these models all exist above a certain critical dimension and the results we see are independent of dimension.
iii.2 L Norm Results
We now look at dynamical transitions in 2D, but using the L norm metric rather than the Euclidean distance. There is an approximately equivalent power law distribution at the critical time, as shown in Fig. 7 for and for varying , in the range: at the same density, . When L norms are used to check for hard sphere overlap, some amount of interpenetration of particles is permitted, as can be seen by the fact that an Lnorm can exaggerate the importance of dominant vector components regardless of the contribution from remaining components. Models of granular media based on L norms have been used to study collisions between nonspherical particles gamba ().
iii.3 The Coalescence Case
Coalescence dynamics are relevant to many physical processes such as nanoparticle and colloidal growth. We consider only dynamical systems in which the total energy is conserved. Our investigation here focuses on whether the trends observed in the Sinai billiard case, namely the behavior of the empirical cluster distribution, extend to the coalescence billiard case.
The empirical cluster mass distribution at has a power law at , as seen in Fig. 8 for the model where all collisions result in coalescence. This is similar to what we observed in the elastic billiards case. For the coalescing cluster distributions, we only considered clusters which have been involved in at least one coalescing collision (of mass ). Billiards which have not coalesced are not considered, as they do not satisfy the neighbor relationship with any other billiards. The coalescing cluster distribution plot compares the magnitude of each cluster [] against its empirical probability.
As increases, also increases due to the decreasing fraction of successful coalescing collisions. For higher , fewer small clusters are formed and large clusters dominate throughout the trial. Figure 9 shows these distributions over multiple runs for different values of . The critical exponent fit is consistent across all models, indicating a large degree of universality in the system. One possible reason for this universality is that all models are at the critical dimension. Such a behavior has been observed by Kang in the case of the Smoluchowski equation and kang () in the elastic case. Our work could be further extended to investigate the form of the kernel for Smoluchowski equation for this model.
The analysis of dynamical phase transitions is relevant to coalescing clusters. Coalescent events have been studied recently in nanoparticle growth trajectories zheng (). This study may be useful in modeling such phenomena.
Iv Conclusion
We have extended the analysis of dynamical phase transitions to higher dimensions, densities, and norms. We have also considered two cases of collisions: elastic collisions and coalescing billiards, and for the elastic case the analysis included the effects of an L norm. We have found universality in the form of a power law for the probability distribution of cluster sizes, with the same critical exponent describing these systems. In the noncoalescing case, the expectation value of the critical time was shown to be determined mainly by the average time between collisions. These observations can be used to predict the onset of criticality and could be used for applications in chemistry such as gas dynamics or polymerization.
References
 [1] Ya.G. Sinai. Construction of dynamics for onedimensional systems of statistical mechanics. Theoret. Math. Phys., 11(2):487–494, 1972.
 [2] Ya. G. Sinai. The construction of cluster dynamics for dynamical systems of statistical mechanics. Moscow Univ. Math. Bull., 29:152–158, 1974.
 [3] N.N. Bogolyubov. Stochastic processes in dynamical systems. Sov. J. Particles Nucl. (Engl. Transl.); (United States), 9, 1978.
 [4] A. Gabrielov, V.I. KeilisBorok, Ya.G. Sinai, and I. Zaliapin. Boltzmann’s Legacy. European Mathematical Society, 2007.
 [5] V.I. KeilisBorok and A.A. Soloviev. Nonlinear dynamics of the lithosphere and earthquake prediction. Springer, 2002.
 [6] A. Gabrielov, V. KeilisBorok, I. Zaliapin, and W.I. Newman. Critical transitions in colliding cascades. Phys. Rev. E, 62(1):237–249, Jul 2000.
 [7] A. Gabrielov, I. Zaliapin, W.I. Newman, and V.I. KeilisBorok. Colliding cacscades model for earthquake prediction. Geophys. J. Int., 143:427–437, May 2000.
 [8] I. Zaliapin, V.I. KeilisBorok, and M.Ghil. A Boolean delay equation model of colliding cascades. part ii: Prediction of critical transitions. J. Stat. Phys., 111(34):839–861, May 2003.
 [9] I. Zaliapin, V.I. KeilisBorok, and M. Ghil. A Boolean delay equation model of colliding cascades. part i: Multiple seismic regims. J. Stat. Phys., 111(34):815–837, May 2003.
 [10] V.I. KeilisBorok, J. H. Stock, A. Soloviev, and P. Mikhalev. Prerecession pattern of six economic indicators in the usa. J. Forecast, 19:65–80, 2000.
 [11] I. Zaliapin, E. FoufoulaGeorgiou, and M. Ghil. Transport on river networks: a dynamical approach. arXiv:0902.1554v1 [physics.geoph], 2009.
 [12] A. Gabrielov, V. KeilisBorok, S. Olsen, and I. Zaliapin. Predictability of extreme events in a branching diffusion model. arXiv:1003.0017v1 [physics.geoph], 2010.
 [13] K. Kang and S. Redner. Fluctuation effects in Smoluchowski reaction kinetics. Phys. Rev. A, 30(5):2833–2836, 1984.
 [14] P. Ghosh. A comparative study of the filmdrainage models for coalescence of drops and bubbles at flat interface. Chem. Eng. Technol., 27(11):1200–1205, 2004.
 [15] M.V. Smoluchowski. Versuch einer mathematischen theorie der koagulationskinetik kollider losungen. Z. Phys. Chem., 92(2):129–168, 1917.
 [16] M.V. Smoluchowski. Drei vortrage uber diffusion, brownsche bewegung und koagulation von kolloidteilchen. Physik. Zeit., 17:557–585, 1916.
 [17] H. Wang, A. Z. Zinchenko, and R. H. Davis. The collision rate of small drops in linear flow fields. J. Fluid Mech., 265:161–188, 1994.
 [18] A. K. Chesters. The modelling of coalescence processes in fluidliquid dispersions: a review of current understanding. Chem. Eng. Res. Des., 69(A4):259–270, 1991.
 [19] G. R. Zeichner and W. R. Schowalter. Effects of hydrodynamic and colloidal forces on coagulation of dispersions. J. Colloid Interface Sci., 71(2):237–253, 1979.
 [20] Ya.G. Sinai. Construction of cluster dynamics for dynamical systems of statistical mechanics. Proc. of Moscow State University, N1, 1974.
 [21] I. M. Gamba, V. Penferov, and C. Villani. On the Boltzmann equation for diffusively excited granular media. Comm. Math. Phys., 246:503–541, 2004.
 [22] A. Hasmy and R. Jullien. Percolation in clustercluster aggregation processes. Phys. Rev. E, 53(2):1789–1795, 1996.
 [23] Ya.G. Sinai and Yu. M. Sukhov. Existence theorem for solutions of the Bogolyubov equations. Theor. Math. Phys, 19(3):560–573, 1974.
 [24] M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato. Packing hyperspheres in high dimensional Euclidean spaces. Phys. Rev. E, 74, 2006.
 [25] H. Zheng, R. K. Smith, Y. Jun, C. Kisielowski, U. Dahmen, and A. P. Alivisatos. Observation of single colloidal platinum nanocrystal growth trajectories. Science, 324(5932):1309–1312, 2009.