Fractal aggregation in charged granular gases
Abstract
The empirical observation of aggregation of small solid particles under the influence of electrostatic forces lies at the origin of the theory of electricity. The growth of clusters formed of small grains underpins a range of phenomena from the early stages of planetesimal formation to aerosol and industrial processes. However, the collective effects of Coulomb forces on the nonequilibrium dynamics and aggregation process in a granular gas, a representative of the above physical systems, are unclear. Here, we establish the premise of a hydrodynamic theory of aggregating granular gases that exchange charges upon collisions and interact via the longranged Coulomb forces. We derive the governing equations for the evolution of granular temperature, charge variance, and number density for a homogeneously aggregating gas. We find that, once the aggregates are formed, the system obeys a physical constraint of nearly constant dimensionless ratio of electrostatic energy to granular temperature ; this constraint on the collective evolution of charged clusters is confirmed both by the theory and the detailed molecular dynamics simulations. The inhomogeneous aggregation among monomers and clusters in the mutual electrostatic field statistically proceeds in a selfsimilar manner. Our work illustrates the defining role of longranged interactions and can lead to novel designing principles in particulate systems.
The electrostatic aggregation of small particles is among the oldest scientific observations. Caused by collisional or frictional interactions among grains, large amounts of positive and negative charges can be generated. These clusters have farreaching consequences: from aerosol formation to nanoparticle stabilization Castellanos (2005); Schwager et al. (2008), planetesimal formation, and the dynamics of the interstellar dust Wesson (1973); Harper et al. (2017); Brilliantov et al. (2015); Blum (2006). The processes accompanying granular collisions and charge buildup can also lead to catastrophic events such as silo failure, or dust explosions.
Experimental investigations of the effects of tribocharging date back to Faraday, and recent in situ investigations have revealed important results Lee et al. (2015); Yoshimatsu et al. (2017); Poppe et al. (2000); Haeberle et al. (2018). However, careful experiments are plagued by a myriad of technical difficulties that often impede unambiguous interpretation Spahn and Sei (2015) of the experiments. A consequence of these difficulties is the lack of consensus about whether electrostatics facilitate or hinder the aggregation process of a large collection of granular particles Spahn and Sei (2015). The aggregation of micronsize particles has been described by meanfield rate equations without considering dissipation or charge transfer mechanisms Ivlev et al. (2002); Dammer and Wolf (2004). Monopolar dilute systems without charge distribution have also been studied Müller and Luding (2011). Intensive studies have been conducted on aggregating systems with shortrange potentials Ulrich et al. (2009); Brilliantov et al. (2018).
Despite its farreaching significance, a statisticomechanical description of aggregation in a granular system with a mechanism of charge transfer is still lacking. The theoretical treatment requires both dissipation of kinetic energy described by a monotonic dependence of the coefficient of restitution on velocities , and also inclusion of longrange electrostatic forces due to the dynamicallychanging charge production. Understanding the growth of charged aggregates requires a statistical approach due to the different kinetic properties and aggregate morphology.
In this work, we derive a Boltzmann equation for the inelastic collisions of grains that interact via Coulomb forces, and produce charges upon collision. We derive the hydrodynamic equations for the number density of the aggregates, the granular temperature , and the charge fluctuations . We find that the dynamics of the charged granular gas approach, but do not overcome, a limiting behavior marked by the value of the dimensionless ratio of electrostatic energy to thermal energy.
To bolster our results, and explicitly consider microscopic fluctuations and morphological structures, we also use threedimensional molecular dynamics (MD) simulations that explicitly include Coulombic interactions and a chargeexchange mechanism. We find that the granular dynamics agree quantitatively with the predictions of the Boltzmann equation. The cooling gas undergoes a transition from a dissipative to an aggregative phase marked by a crossover in the advective transport.
.1 Kinetics
To obtain a general description, we consider the single particle chargevelocity phase space probability density function , where the particle velocity , charge , and position are the phase space variables, and denotes time. We consider a homogeneously aggregating (homogeneous and quasimonodisperse) granular gas with charge exchange and dissipation upon collisions. In this limit, the spatial and particlesize dependence of drops out and its time evolution is given by the simplified Boltzmann equation Pitaevskii and Lifshitz (2012); Brilliantov and Pöschel (2010)
(1) 
where we define as the generalized collision integral which includes dissipation as well as charge exchange during particle collisions. We will now elucidate how the charge exchange mechanism modifies the collisions.
Let us consider contact collisions of particles and with precollision velocitycharge values and , respectively. In the ensemble picture, particle collisions will change in the infinitesimal phasespace volumes and , centered around and , respectively. The number of direct collisions per unit volume which lead to loss of particles from the intervals and in time are
(2) 
where , is the unit vector at collision pointing from the center of particle towards particle , and is the differential collision crosssection. The Heaviside step function selects particles coming towards , while we use to ensure that a contact with an approaching particle takes place only when the Coulomb energy barrier can be overcome, where , F m is the vacuum permittivity, and is the particle diameter. If the interaction is repulsive, is positive, and only if . In case of attractive interaction, is negative and thus always. Essentially, filters repulsive interactions which do not lead to a physical contact between particles.
Consider now particles with initial velocitycharge values and in the intervals . The number of particles per unit volume which, postcollision, enter the interval and in time are
(3) 
The net change of number of particles in time per unit volume, then reads (see Methods for details)
(4) 
where and are Jacobians of the transformations for and respectively (see Methods), which lump together the microscopic details of the collision process, namely dissipation and charge exchange in the present study.
Integrating over all incoming particle velocities and charges from all directions, dividing by and taking the limit , we obtain the formal expression for the collision integral
Here we assume that the differential collision cross section, and the contact condition specified by , retain their form for direct and inverse collisions. At this point the particle encounters which do not lead to a physical contact have been excluded using . While taking moments of , a fraction of those contact collisions that lead to aggregation is accounted for by taking the limit for certain conditions on the relative velocity , and by considering the charge transferred to particle equal to the charge on particle [Methods, Eq. (24)(29)]. In , distant encounters, which do not lead to a contact between particles (glancing collisions) are neglected and the charge exchange and dissipation is considered only during the contact. The longrange effect is incorporated via the collision crosssection [see Eq. (34)].
After setting up the collision integral, we derive the macroscopic changes of number density , granular temperature , and the charge variance (due to aggregation, dissipation and charge exchange, respectively) by taking the moments of the Boltzmann equation (see Methods: splitting restitution and aggregation). The particles are initially neutral and the charge on them is altered either by collisions or during aggregation. However, due to charge conservation during collisions and aggregation, the system remains globally neutral and the mean charge variation is zero. An alternative choice is thus . In order to obtain closed form equations, and for analytical tractability, we assume quasimonodispersity and homogeneity of the aggregating granular gas at any given time, as illustrated in Fig. 1. This means that during aggregation the mass of the clusters is assumed to grow homogeneously, while their numbers decreases in a given volume. In the following we show that even under these assumptions our hydrodynamic equations provide predictions close to the MD results. We also show that the time evolution of , however, is more involved.
After integration under the assumption of uncorrelated, Gaussian charge and velocity distributions (see Methods) we find the governing equations
(5) 
(6) 
(7) 
where the coefficients , , , , (listed in Table 1) are functions of the time dependent ratio
(8) 
between electrostatic energy and granular temperature. We term the ratio as Bjerrum number. The variables are also timedependent functions of aggregating cluster mass , size , and the material constants (Table 1).
First of all, we note the consistency of the above set of equations with the classical Haff’s law in the absence of collisional charging. In this limit , , and we obtain , and , reducing to Haff’s law.
For nonzero charge exchange, the coupled equations (5)(7) are solved numerically. Figure 2 shows the temporal evolution of the hydrodynamic quantities , , , as well as the dimensionless parameter . Initially, the collisional velocities remain larger than the threshold (see Eq. (23) and (39) in Methods for the role of the threshold ). In this time regime, the collisions are purely restitutive leading to charge exchange and dissipation. The charge exchange increases while the dissipation reduces . The number density , and thus size of the effective particles, remain unchanged as there is no aggregation. The Bjerrum number , however, increases. As approaches the threshold , , and the aggregation initiates, the concentration begins to decrease and thus the size of the aggregates (now considered as new “‘effective” particles; see Fig. 1) increases. The individual particles, or monomers, cluster in such a way that the charge variance on the cluster population now begins to reduce. The theory shows that once aggregation sets in, the aggregating granular gas obeys the constraint
(9) 
This physical limit is also confirmed by granular MD simulations (see Methods, and Singh and Mazza (2018), for details on MD). At long times, adjusts itself such that the right hand sides of the equations for number density, temperature and charge variance [Eq. (5)(7)] remain real valued.
To understand the origin of this approximate conservation law in the theory, we solve Eq. (5)(7) for a system with only dissipative collisions and without aggregation; hence, the cluster size remains unchanged. The results are shown in Fig. 2 (dashed lines). In this limit of only restitutive kinetics, increases continuously, while deviates from Haff’s law. If the complete aggregation dynamics are considered (i.e., dissipation and aggregation) we find that Haff’s law for aggregating clusters is closely obeyed. This is also confirmed by our MD simulations.
Figure 2(b) shows that, after the initiation of aggregation, the charge variance of the aggregating clusters also follows a power law decay. It is notable that a different charge exchange model might provide a different charge buildup rate during the purely dissipative (restitutive) phase, however, the decay of charge variance during the aggregation of the monomers into the clusters is not expected to be influenced by a change in the charge exchange mechanism. The reason is that aggregation sets in at relatively low temperature where the relative thermal motion of monomers inside the clusters is insignificant, and thus collisional chargeexchange is expected to be negligible. In our theory, the charge transferred to a particle during aggregation with particle is equal to the charge on particle as they are considered merged into one single particle (see Fig. 1). Thus the charge variance of cluster population is reduced by the aggregation process, rather than by the collisional charge exchange in the aggregation phase. In the MD simulations, the decay of during the aggregation phase agrees well with the theory, even though we observe that the charge exchange in the simulations leads to a symmetric but nonGaussian charge distribution during the restitution phase (see Methods, Fig. 5). Although neglected in the current model, the discharging of the monomers inside the aggregates at very low collision speeds is neglected, which might change the charge variance decay rate.
The key finding here is that the decay of the charge variance and the growth of the size of the aggregates is balanced by the decay of temperature during the aggregation, resulting in the stationary value of . The steady state of is apparent in the theory. The granular MD simulations suggest and confirm the upper limit of . We simulate the above results under different initial filling fractions and find that the constraint persists with variation of . The temperature of the cluster population also closely follows Haff’s law for different (see Fig. 2 insets).
The number density evolution from the MD simulations is more involved than the prediction of the quasimonodisperse kinetic theory. The theory indicates the decay of cluster density only in an average sense (Fig. 2(c) and insets). To explore the mechanism behind the nonmonotonic evolution of , we first explore the cluster dynamics and nature of the structures from the full spatiallyheterogeneous MD simulations.
.2 Inhomogeneous aggregation and fractal growth
To gain access to the spatial structure formation in the gas, we perform a detailed cluster analysis of the results from granular MD simulation (see Methods and Singh and Mazza (2018) for more details); our cluster analysis code is made publicly available (see Support Information).
Figure 3 shows snapshots of the clusters at different times from a typical run.
The morphology of the aggregates is studied by computing the average fractal dimension Mandelbrot (1977); Jullien (1987) for the cluster population from the scaling law between cluster mass and radius of gyration , where the index runs over all monomers in a given aggregate. Figure 4(ac) shows scatter plots for versus at different times. Figure 4(d) shows the time evolution of the exponent for varying filling fraction. The average fractal dimension lies in between the average values reported for the ballistic clustercluster aggregation (BCCA, ) and the diffusionlimited particlecluster aggregation (DLPCA, ) Blum (2006); Smirnov (1990) models. In time, is dynamic and changes across the two model limits. These results indicate that the aggregate structures retain their fractal nature over time.
The BCCA and DLPCA mechanisms of aggregation have been reported in neutral dust agglomeration (e.g hit and stick, ballistic motion, Blum (2006)), wet granulate aggregation (sticking due to capillary bridges and ballistic motion, Ulrich et al. (2009)), colloidal aggregation (van der Waals and repulsions Lebovka (2012)), and hit and stick agglomeration in Brownian particles under frictional drag Kempf et al. (1999). The observation that lies between the reported average values of for BCCA and DLPCA indicates the presence of mixed characteristics from both of these simplified models. The size distribution in an aggregating charged granular gas (see Singh and Mazza (2018)) tends to resemble a DLPCAlike behavior where the smaller size aggregates are larger in number, in contrast to a BCCAlike model where the size distribution is typically bellshaped Blum (2006). On the other hand, the monomer motion is found to be highly subdiffusive Singh and Mazza (2018) in agreement with the BCCA model. In addition, the Coulombic interactions will cause considerable deviations from the shortranged or ballistic propagation typical of the BCCA or DLPCA models. The longrange forces due to a bipolar charge distribution lead to the intermediate value of between what is observed in the above two physical aggregation models.
.3 Interplay between fractals and macroscopic flow
Apart from the longrange effects, the morphology of the aggregates is also altered by additional mechanisms. We discuss two physical processes that are not captured in the analytical theory, but that we investigate via our MD simulations.
First, in our Boltzmann kinetic description, the collisions between aggregates at any given time are considered as collisions between two spheres with sizes equal to the average size of the aggregate population. This is a quasimonodisperse assumption typically used in clustercluster aggregation (CCA) models. The quasimonodispersity however neglects the morphology and surface roughness of the colliding aggregates. Collisions between aggregates with large size difference, and with individual monomers, are also simplified.
Secondly, granular gases are characterized by the emergence of a convective flow Hummel et al. (2016) which we find, in the present case, inducing the nonmonotonicity in the temporal evolution of the number density [Fig. 2(c)]. To quantify the emergence of the macroscopic flow, we calculate the Mach number
(10) 
Here the system is divided into equal sized cubic boxes and is the center of mass velocity in a box (or local advective velocity), scaled by , the thermal velocity of the monomers in the gas. The time evolution of is shown in Figure 2(c) (inset), which indicates that grows at a higher rates when the number density evolution becomes nonmonotonic (), indicating a coupling between aggregation and the generation of macroscopic flow. Due to the macroscopic flow (see also the Supplementary Media), aggregates which are weakly connected are prone to fragmentation. This results in an intermediate regime where the concentration increases instead of decreasing.
Excluding the two mechanisms above explains the shortcoming of our quasimonodisperse Boltzmann theory in capturing the nonmonotonic behavior of occurring at , and only exhibits the general downward trend of the number density [Fig 2(c)].
.4 Conclusions
We have derived the rate equations for the evolution of the number density , granular temperature , and charge variance of the cluster population in a charged, aggregating granular gas. In contrast to wellknown Smoluchowskitype equations, we have explicitly coupled to the decay of and charge variance. We have compared the results with threedimensional molecular dynamics simulations and the outcomes of a detailed cluster analysis, and have explored the morphology of the aggregating structures via fractal dimension.
Taken together, our results indicate that the aggregation process in a charged granular gas is quite dynamic, while respecting some physical constraints. The growth process obeys , while morphologically, the clusters exhibit statistical selfsimilarity, persistent over time during the growth. The fractal dimension and growth of structures is intermediate between the BCCA and DLPCA models. We also demonstrate that the application of a purely dissipative kinetic treatment is not sufficient to make predictions about global observables such as and in an aggregating charged granular gas.
Finally, we believe that our kinetic approach can be applied to study aggregation processes in systems such as wet granulate collections with ion transfer mechanism Lee et al. (2018); Zhang et al. (2015), dissipative cell or active particle collections under longrange hydrodynamic and electrostatic effects Yan et al. (2016); Friedl and Gilmour (2009), and charged iceice collisions Dash et al. (2001).
.5 Acknowledgements
We gratefully acknowledge highly constructive comments by Jürgen Blum, Eberhard Bodenschatz, Stephan Herminghaus, and Reiner Kree. We thank the Max Planck Society for funding.
.6 Methods and supporting information
Kinetics and generalized collision integral. Let us consider the number of particles per unit volume having initial velocitycharge values and in the intervals and which, postcollision, enter the in the intervals and in time are
(11) 
and thus the net increase of number of particles per unit time and volume is . We can relate the primed variables to the unprimed via
(12) 
where is the Jacobian of the transformation for viscoelastic particles. Here is a material constant. To obtain the transformation , we consider the ratio of relative charges after and before the collision
(13) 
and in addition we impose charge conservation during collisions
(14) 
The above two relations finally provide the transformation
(15) 
where, for example, for a constant . This means that the differential chargespace volume element shrinks or expands by a factor of . In general, for velocity and chargedependent case the expressions of and can be quite complicated as it depends on how the charge exchange takes place during collisions and its dependence on myriad factors (such as size, composition, and crystalline properties). Incorporating the above phasespace volume transformations due to collisions, the net change of number of particles per unit phasespace volume and in time reads
(16) 
where we assume that the differential crosssection and the contact condition specified by are the same for direct and inverse collisions. Finally, dividing by , and integrating over all incoming particle velocities and charges from all directions in the limit , we obtain a formal expression for the collision integral
(17)  
At this point the particle encounters which do not lead to a physical contact have been excluded using , however, collisions that lead to aggregation have not been explicitly accounted. We do so by taking the limit for certain conditions on the relative velocity , and by considering the charge transferred to particle equal to the charge on particle ((24)(29)). In , distant encounters, which do not lead to a contact between particles (glancing collisions) are neglected and the charge exchange and dissipation is considered only during the contact. The longrange effect is incorporated via collision crosssection in Eq. (34).
Splitting restitution and aggregation.
The time rate of change of the average of a microscopic quantity is obtained by multiplying the Boltzmann equation for by and integrating over , i.e.
(18) 
It can be shown that
(19) 
where and is the change of during the collision between pair , and the prime denotes a post collision value. We note that the transformations in Eq. (12) and (15) are reversed back while integrating weighted with quantity of interest . We consider the number density, the kinetic energy or granular temperature, and the charge variance (the system is globally neutral and the mean charge variation is zero), respectively
(20)  
(21)  
(22) 
At this point we differentiate the restitutive or dissipative collisions from aggregative ones by splitting as
(23) 
where . If , the particles collide and separate after the collision irrespective of the sign of (attractive or repulsive). This leads to dissipation of energy with finite nonzero , and charge exchange according to a specified rule. The aggregative part is zero in this case. If and (attractive encounters at low velocities), the particles collide and aggregate with , and with charge exchange to particle equal to . If and (repulsive encounters at low velocities), no physical contact takes place between the particles which leads to neither dissipation nor aggregation ().
The expressions for and are obtained as follows. The particle number does not change during a dissipative collision but reduces by one in an aggregative collision, i.e.
(24) 
For the granular temperature,
(25) 
where we take the limit for the aggregation.
The change in the charge variance is obtained as
(26) 
where equals the charge transferred to particle during its collision with particle , and is the mean charge on the pair. For the charge transfer, we consider
(27) 
which together with Eq. (26), gives
(28) 
For aggregation, the charge tranferred to particle equals the charge on the merging particle , i.e, , which gives
(29) 
Putting Eq. (24),(25),(28), and (29) in (23), and then (23) in (19), the resulting integrals are solved (see Appendix for details), assuming the statistical independence of chargevelocity distribution function, i.e., , and assuming that they are Gaussian. In addition to the charge exchange, the coefficient of restitution is taken as velocity dependent, i.e
(30) 
while the longrange effects due to Coulomb interactions are taken into account by the change in collision cross section.
Granular MD simulations. The equation of motion of the form
(31) 
is solved for a setup with periodic boundary conditions in a cubic box of size , where is monomer diameter. Here the elastic force strength is , and a dissipative force strength of is set to obtain clustering in a finite size system (See Singh and Mazza (2018) for more details). Also . The Coulomb force strength is varied between . The longrange Coulomb forces for a setup with periodic boundary conditions is challenging and conditionally convergent as it depends on the order of summation. We employ the Ewald summation that converges rapidly, and has a computational complexity . The algorithm is partially parallelized and highly optimized on graphics processing unit (GPU). In our simulations, the total computing time to reach simulation time for a typical simulation with , including the longrange electrostatic forces, is of the order of weeks. See Singh and Mazza (2018) for more details.
Codes and media. The following computer codes and media are made public:

MATHEMATICA framework to solve the restitutive and aggregative parts of the kinetic integrals.

The cluster analysis code in MATLAB to obtain the fractal dimension, cluster size distribution, average cluster size, and other statistical quantities, is made available at https://gitlab.com/cphyme/clusteranalysis.

Movie of the evolving and aggregating structures in the charged granular gas.
Data availability. All data are available from the authors.
Dimensionless quantities to laboratory values. Typically nonBrownian growth in planetary dust is dominant for monomer sizes of several Zsom et al. (2010) and the growth barrier problem Spahn and Sei (2015) starts to arise near . The mass of silica particles in this range of size is . Thus setting gives a velocity reference value . This fixes the time reference . Thus in our results the growth over units of nondimensional time approximately implies growth over , or for such size and mass particles. The average size in the growth period grows approximately by one order of magnitude (e.g. the growth of the largest cluster is from to in time under these reference values, and for initial filling fraction ).
.7 Restitutive and aggregative parts of (19).
To solve the integrals in Eq. (19) for different from (20)(22), we assume that the normalized velocity as well as charge distribution in the gas at any time remain Gaussian, and the two are uncorrelated, i.e
(32) 
In Figure 5, we show charge distribution of monomers obtained from typical simulation runs. In (32), we assume that shape of the scaled charge distribution remains close to Gaussian during size growth.
The coefficient of restitution is velocity dependent, i.e
(33) 
The attractive or repulsive longrange effects are emulated through an effective differential crosssection for a binary collision, which changes depending upon the sign and magnitude of charges on the particle pair and their relative velocity, according to
(34) 
where is the differential crosssection per unit solid angle . The expression is independent of , and thus the total cross section is . For neutral particles, , and thus . For (repulsive encounters), , while for (attractive encounters), .
Below we explain the solution procedure for the restitutive part of the equation for . Similar procedure can then followed for the aggregative part of , except the fact that in this case and integration limits are selected by , where the shorthand . Plugging Eq. (25) into (23) and then Eq. (23) into (19), we find