Electrification in granular gases leads to constrained fractal growth.

The empirical observation of aggregation of dielectric 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 aerosols. However, the collective effects of Coulomb forces on the nonequilibrium dynamics and aggregation process in a granular gas - a model representative of the above physical processes - have so far evaded theoretical scrutiny. Here, we establish a hydrodynamic description of aggregating granular gases that exchange charges upon collisions and interact via the long-ranged Coulomb forces. We analytically derive the governing equations for the evolution of granular temperature, charge variance, and number density for homogeneous and quasi-monodisperse aggregation. We find that, once the aggregates are formed, the granular temperature of the cluster population, the charge variance of the cluster population and the number density of the cluster population evolve in such a way that their non-dimensional combination obeys a physical constraint of nearly constant dimensionless ratio of characteristic electrostatic to kinetic energy. This constraint on the collective evolution of charged clusters is confirmed both by  our theory and our detailed molecular dynamics simulations. The inhomogeneous aggregation of monomers and clusters in their mutual electrostatic field proceeds in a fractal manner. Our theoretical framework is extendable to more precise charge exchange mechanisms, a current focus of extensive experimentation. Furthermore, it illustrates the collective role of long-ranged interactions in dissipative gases and can lead to novel designing principles in particulate systems.

Introduction

The electrostatic aggregation of small particles is ubiquitous in nature and ranks 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 far-reaching consequences: from aerosol formation to nanoparticle stabilization[1,2], planetesimal formation, and the dynamics of the interstellar dust[3-6]. The processes accompanying granular collisions, charge buildup and subsequent charge separation 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[7-11]. However, technical difficulties plague even careful experiments and often impede their unambiguous interpretation[12]. A source of these difficulties is the lack of consensus about whether electrostatics facilitate or hinder the aggregation process of a large collection of granular particles[12]. Despite considerable effort[13-24] a statistico-mechanical description of aggregation in a dissipative granular system with a mechanism of charge transfer is still lacking. The theoretical treatment requires reconsideration of the dissipation of kinetic energy conventionally described by a monotonic dependence of the coefficient of restitution on collision velocities , and also inclusion of long-range electrostatic forces due to the dynamically-changing 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 present a modified Boltzmann description for the inelastic and aggregative collisions of grains that interact via Coulomb forces, and exchange charges upon collision. We derive the hydrodynamic equations for the number density n, the granular temperature T, and the charge fluctuations of the aggregates under the assumptions of homogeneous and quasi-monodisperse aggregation. We find that the dimensionless ratio of the characteristic electrostatic energy to the granular temperature T, i.e. , approaches but stays below  unity. Here is the Coulomb’s constant while d is the characteristic size of the aggregates.

To bolster our results, and explicitly consider fluctuations in dynamics and morphological structures, we also use three-dimensional molecular dynamics (MD) simulations that explicitly include Coulombic interactions and a charge-exchange 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. We explore the morphological dynamics of the inhomogeneous aggregation via the mean fractal dimension and their interplay with the mesoscopic flow.

Kinetics

In general, agglomeration in a three-dimensional collisionally charging cooling granular gas is a spatially inhomogeneous process which involves the interplay between dissipation, time-varying size distribution of aggregates, charge fluctuations and exchange mechanism during collisions, long-range forces, and collective effects[25]. This complexity is illustrated in Fig. 1 which shows snapshots of cooling clusters from a typical MD simulation, beginning from a homogeneous and neutral state (see Methods and[25] for MD). In the following we establish a modified Boltzmann approach for this intricate dynamics of the aggregation process, which predicts novel physical limits.

Figure 1

Snapshots of the aggregating charged granular gas at different non-dimensional times: (a) , (b) , and (c) . Here clusters containing 10 particles or more () are shown, and each color represents a different cluster. Clusters are identified on the basis of the monomer distances: if the centers of two particles are separated by a particle diameter, or less, they belong to the same cluster. See Methods for non-dimensional time , Coulomb force strength , and other reference scales.

We consider the single particle probability distribution function , where the particle velocity v, charge q, position r and particle size d, are the phase space variables, and t denotes time. We specialize to a homogeneous and quasi-monodisperse aggregation scenario (i.e. the size is assumed to vary in time but spatially mono-dispersed, see schematic representation in Fig. 2). Under these limits, the spatial and particle-size dependence of f drops out, i.e. only, and its time evolution is given by the simplified Boltzmann equation[26,27]valid at any time instant t. Here we define Icoll as the modified collision integral which includes dissipation as well as charge exchange during particle collisions. We will now elucidate how the charge exchange mechanism and particle size growth modify the collisions.

Figure 2

(a) When charged particles collide (particle group I), two possible pathways are considered in the kinetic theory depending on their velocities and charges: (i) a typical inelastic collision (dissipation with velocity dependent coefficient of restitution) with charge exchange (particle group II) if the relative collision speed is above a threshold speed , and (ii) a collision that leads to the aggregation (with coefficient of restitution equal to zero) and merging of charges (particle group III), if the relative collision speed is below and the particles are oppositely charged. The latter event introduces a size distribution of aggregates, which we further simplify in our kinetic theory by adjusting the size of all the particles (particle group IV). Thus, the particle size is assumed to remain monodispersed during the course of aggregation. Primed variables represent post-collision or post-aggregation values. (b) Coefficient of restitution in the present kinetic theory. Notice that the threshold also evolves with time due to charge-exchange and aggregation events, and the coefficient of restitution is zero only if . The above kinetic approach is compared with the MD simulations where the size distribution of aggregates is polydisperse.

Let us consider contact collisions of particles i and j with pre-collision velocity-charge values and , respectively. In the ensemble picture, particle collisions will change in the infinitesimal phase-space volumes and , centered around and , respectively. The number of direct collisions per unit spatial volume which lead to loss of particles from the intervals and in time Δt arewhere , n is the unit vector at collision pointing from the center of particle i towards particle j, and dσ is the differential collision cross-section [Eq. (38)]. The Heaviside step function selects particles coming towards i, while we use to ensure that a contact with an approaching particle takes place only when the Coulomb energy barrier can be overcome, where , is the vacuum permittivity, and d is the particle diameter at time t. 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 velocity-charge values and in the intervals , . The number of particles  (inverse collisions) per unit volume which, post-collision, enter the interval and in time Δt is

The net change of number of particles in time Δt per unit volume, then reads (see Methods)where and are the Jacobians of the transformations for and , respectively, 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 Δt and taking the limit , we obtain the formal expression for the modified 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. The particle encounters which do not lead to a physical contact have been excluded using . While taking moments of Icoll, a fraction of those contact collisions that lead to aggregation is accounted for by taking the limit for certain conditions on the relative velocity v, and by considering the charge transferred to particle i equal to the charge on particle j. In Icoll, 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 long-range effect is incorporated via the collision cross-section [Eq. (38)].

After setting up the collision integral, we derive the macroscopic changes of number density , granular temperature , and the charge variance for a homogeneously aggregating granular gas by taking the moments of the Boltzmann equation (see Methods), in the absence of any macroscopic velocity V. Here Nagg  is  the number of aggregates in a spatial volume V; , and q are the aggregate mass, velocity and charge respectively, while subscript “mean” signifies the mean value of the quantity. 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. The next choice is thus . In order to obtain closed form equations, and for analytical tractability, we assume quasi-monodispersity and homogeneity of the aggregating granular gas at any given time, as illustrated in Fig. 2. This means that during aggregation the mass of the clusters is assumed to grow homogeneously, while their numbers decrease in a given volume.

We assume that the charge and velocity distributions are uncorrelated, and their properly scaled form remains Gaussian. After integration we find the governing equationswhich are coupled via a time-dependent dimensionless ratiobetween charge variance, granular temperature, and aggregate size. The terms g are time-dependent functions of and material constants , (Methods Table 1 and Eq. (32–36)). We term the ratio as Bjerrum number. In Eq. (9), d represents the size of a particle, also evolving with time during aggregation [Fig. 2]. Notice that as f is considered independent of d during aggregation, an explicit equation for d is required. For this we consider the total mass M, system volume V, and particle material density to be constant, which fixes the relation between particle size d and particle number density n, according toand closes the Eqs set (6–8). Here D is the fractal dimension. In the solution of the above analytical equations, we have assumed (spherical aggregates). Below, we will compute values of D using the MD simulations. The above set of equations is consistent with modified Haff’s law for a velocity dependent coefficient of restitution[28-30] in the absence of collisional charging. In this limit , , and we obtain , and , whose solution is the modified Haff’s law.

Table 1

Expressions of the coefficients in Eqs (5–7) in the main text.

Coefficient	Expression
ℬ	k_e〈δq^2〉/(Td)
l	4−ℬ^2
l_1	1−ℬ^2
CresT	802^1/5d^2C_ϵm/(21πm^8/5)
CaggT	d^2/(8πm)
Cresq	1.4080d^2CΔq2(2+η)(4/m)^η+1/2/π
Caggq	4d^2/πm
Caggn	d^2/(2πm)
a1T	ℬ(8−5ℬ^2)l+π(16−10ℬ^2+3ℬ^4)
a2T	−32+20ℬ^2−6ℬ^4
a3T	l(2π(ℬ^2−4)^2+ℬ(−32+28ℬ^2+ℬ^4))
a4T	−4(32−20ℬ^2+9ℬ^4)
a5T	lℬ(−32+28ℬ^2+ℬ^4)
	+2π(ℬ^2(20−8l)+ℬ^4(−9+l)+16(−2+l))
a6T	4(32−20ℬ^2+9ℬ^4)
a1q	2ℬl(ℬ^2−1)−π(4−2ℬ^2+ℬ^4)
a2q	2(4−2ℬ^2+ℬ^4)
a3q	16lπl_1(4−5ℬ^2+ℬ^4)^2
a4q	−πℬl(π(ℬ^2−4)(1+2ℬ^2)
	+2ℬl_1(−8−54ℬ^2+33ℬ^4+2ℬ^6))
a5q	2πℬl^5(1+2ℬ^2)
a6q	16πℬl15(4+5ℬ^2)
a7q	8lπl_1(4−5ℬ^2+ℬ^4)^2
a8q	−πℬl^5(1+2ℬ^2)
a9q	πℬl_1ℬl(8+54ℬ^2−33ℬ^4−2ℬ^6)
a10q	8πℬl_1(ℬ^2−1)^2(4+5ℬ^2)
a1n	−ℬ^2l^3
a2n	l_1lℬ(−4+7ℬ^2)+l_1π(ℬ^2−1)(8−4l+ℬ^2(−6+l))
a3n	4l_1(4−7ℬ^2+3ℬ^4)
a4n	−ℬ^2l^3
a5n	l_1l(ℬ(−4+7ℬ^2)+π(4−5ℬ^2+ℬ^4))
a6n	−4l_1(4−7ℬ^2+3ℬ^4)

Here and are the mass and size of the aggregates. The material constant [Eq. (31)] influences the viscoelastic properties of the particles, while and [Eq. (28)] influence the charge buildup. Other notations are as described in the main text.

Results

The time evolution of T, and n of aggregate population from MD simulations is shown in Fig. 3. The aggregate temperature from simulations is  computed as . Notice that v and m are center of mass velocity, and mass of the k aggregate respectively, and should not be confused with monomer velocities and masses. V is the local advective velocity in the neighborhood of k aggregate. Similarly, , and , where Nagg is the total number of aggregates and V is the system volume. In the MD simulations, oppositely charged monomers hold together into mechanical contact when the elastic repulsion between them balances the Coulomb attraction. This is the scenario at low granular temperature when the dissipative force term, and the inertial term, in the equation of motion [Eq. (56)] is negligible compared to the elastic and Coulomb terms. We identify the aggregates with the following condition: if the distance between the centers of two monomers r is less than or equal to the monomer diameter d0, they belong to the same aggregate. Once all the aggregates in the system are identified using the above definition, the aggregate velocity and mass as a whole is computed. If v are the velocities of monomers in k aggregate, the aggregate velocity v is computed by , which is the center of mass velocity of the aggregate. The mass of the aggregate is simply . Similarly, the net charge on k aggregate is .

Figure 3

Evolution of (a) temperature of cluster population, (b) charge variance of cluster population, and (c) number density of cluster population, for different monomer filling fractions and charge strength . (d) The granular temperature, charge variance and average size of the cluster population during aggregation evolve in such a manner that their non-dimensional combination (see also Fig. 4). Both temperature and charge variance of cluster population decay as power laws (with exponents marked as legends). The number density evolution, however, is highly dynamic and exhibits a non-monotonic behavior due to emergence of mesoscopic flow (see Tables 2 and 3 for , , , ).

Figure 4

The granular temperature, charge variance and average size of the cluster population during aggregation evolve in such a manner that their non-dimensional combination . This is not captured in the kinetic theory if only restitutive (no aggregation) collisions are considered (dashed line). The granular MD simulations (symbols) confirm the analytical results. (inset for different monomer filling fraction ).

Table 2

Simulation parameters.

Parameter	Expression	Value
N_mon	No. of monomers	8516, 19652, 50016, 89746
L^3	System size	70×70×70
φ	Monomer filling fraction	0.013, 0.030, 0.076, 0.137
ℰ	αlref3/2t_refm_refv_ref	278.0
D	βlref1/2t_refm_ref	27.8
K	k_eqref2t_refm_refv_reflref2	0.4, 1.0, 5.0

Table 3

Reference values in SI units for conversion of non-dimensional results to laboratory relevant values.

Reference scale	Description	Value
m_ref	Particle mass reference	1.52×10^−4kg
l_ref=d_0	Length reference=monomer diameter	4.78×10^−3m
v_ref	Velocity reference	1.0 ms−1
t_ref=d_0/v_ref	Time reference reference	4.78×10^−3s
T_ref=13m_refvref2	Temperature reference	0.51×10^−4kgm^2s^−2
k e	Coulomb’s constant	8.98×10^9Nm^2C^−2
α	Elastic constant of particles	2.67×10^4kgm^−1/2s^−2
β	Viscous constant of particles	1.28×10^1kgm^−1/2s^−1
q ref	Charge reference	2.0×10^−8C–5.6×10^−9C

Initially (), the relative monomer collision velocities v remain larger than the time varying threshold (see Methods for the treatment of threshold b in the kinetic theory). In this time regime, the collisions are primarily restitutive, leading to either Coulomb scattering without collision, or charge exchange and dissipation without considerable aggregation. The dissipation reduces T [Fig. 3(a)] while the charge exchange increases [Fig. 3(b)]. In this time regime the number density n, and thus size d of the aggregates, is altered only moderately due to those low velocity attractive monomer encounters which lead to aggregation [Fig. 3(c)].

As a result of our kinetic formulation, the dynamics of n, and can be collected into the evolution of , shown in Fig. 3(d). The Bjerrum number initially increases, which indicates that temperature decreases at a faster rate than the rate of increase of charge variance and the aggregate size. As the relative velocities v approach the threshold , . Near this time, the dynamics cross over to aggregative collapse. The individual particles, or monomers, cluster in such a way that the charge variance of the cluster population now begins to reduce. The temperature of the aggregate population keeps decreasing at the same rate with a slight dip near the aggregative collapse. The number density starts to evolve non-monotonically. We explore the non-monotonicity in the next section. These results are robust under variation of the initial monomer filling fraction and the charge strength [Fig. 3].

After the initial time regime and the aggregative collapse, the crossover in the dynamics is depicted in Fig. 4, where the evolution of , and its comparison with the solution of Eqs (6–8) is highlighted. We solve the full kinetic equations Eqs (6–8) including the aggregation kinetics [Fig. 4 (solid line)]. We also solve Eqs (6–8) for a system with only dissipative collisions and without aggregation; hence, the cluster size d remains unchanged. These results are shown in Fig. 4 (dashed line). In this limit of only restitutive kinetics, increases continuously above the limiting value 1. The purely restitutive kinetic theory thus fails to predict the MD results. When aggregation is explicitly treated (solid line), the theory predicts an upper limit during the growth. The theory shows that once aggregation sets in, the aggregating granular gas obeys the constraint

The upper physical limit predicted in the theory, , is endorsed by the granular MD simulations under moderate variation of and . It is notable that at later times, the limit allows the right hand sides of the equations for number density, temperature and charge variance [Eqs (6–8)] to remain real valued during the aggregation process. This mathematical indication confirms the effectiveness of the quasi-monodisperse picture [Fig. 2] considered in the present study.

The initiation of aggregation brings about a power law decay in the charge variance [Fig. 3(b)]. 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 aggregative phase is not expected to be influenced by charge exchange mechanism. The reason is that aggregation sets in at relatively low temperature where the motion of monomers, if any, inside the clusters is significantly decreased, and thus collisional charge-exchange is expected to be negligible. If two oppositely charged particles i and j collide with a speed which is below the aggregation threshold speed, it is considered that they form a single aggregate with net charge on it [Fig. 2]. Thus, the charge variance of the cluster population is reduced by the aggregation process, rather than by the collisional charge exchange.

Our key theoretical finding is that after the aggregative collapse, the decay of the charge variance of aggregate population 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 constraint is robust in the theory, while the granular MD simulations suggest and confirm the upper limit of . It is also intriguing that the temperature of the cluster population still closely follows modified Haff’s law for different and , despite complex heterogeneous aggregation-fragmentation events and the long-range electrostatic interactions.

The number density’s temporal evolution obtained from the MD simulations reveals a more intricate non-monotonic dynamics.  Initially  it decreases during small aggregate formations due to low velocity attractive monomer encounters. In an intermediate time regime, the emergence of mesoscopic particle fluxes triggers fragmentation events and the aggregate numbers increase. We quantify the emergence of mesoscopic flow using the Mach number (see Methods and Media therein). After this intermediate time regime, the aggregation again takes over and the number density of clusters starts to decrease. The non-monotonic evolution of n causes a dip in after the aggregative collapse () [Figs 3 and 4]. For temperature values in this regime, the charge transfer events between monomers are statistically ineffective, and the evolution of is primarily dominated by aggregation events. We find that a maximum of charge fluctuations occurs near this crossover.

The size difference between aggregates in the MD simulations further adds to the complexity of ’s evolution after the aggregative collapse, which is neglected in the homogeneous and quasi-monodisperse aggregation kinetic theory. However, the theory still clearly predicts the growth of and selects a unique upper limit after the aggregative collapse. To further explore the mechanisms behind the non-monotonic evolution of n, we explore the spatially heterogeneous cluster dynamics and nature of the structures from the MD simulations.

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 Fig. 1.

The morphology of the aggregates is studied by computing the average fractal dimension [31,32] of cluster population from the scaling relation between cluster masses m, and radii of gyration , where the index i runs over total number of monomers Nmon in a given aggregate. Here r is the position vector  of the i monomer in the given aggregate and rmean is the position vector  of the center of mass of that aggregate. Once m and R of all the aggregates in the system are calculated, we compute the R versus m scatter-plots at different times, for example in Fig. 5(a–c). We repeat this for different times, and calculate as the slope of the fit to versus scatter-plots. Figure 5(d) shows the time evolution of the exponent for  different filling fraction. The average fractal dimension lies between the average values reported for the ballistic cluster-cluster aggregation (BCCA, ) and the diffusion-limited particle-cluster aggregation (DLPCA, )[6,33] 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.

Figure 5

(a–c) The scaling between cluster mass and their radius of gyration , at different times, and (d) the time evolution of , thus obtained, for different filling fractions. The average fractal dimension in the aggregating charged gas varies  between the average values reported for ballistic cluster-cluster aggregation (BCCA, ) and diffusion-limited particle-cluster aggregation (DLPCA, )[6,33]. .

The BCCA and DLPCA are popular models for aggregation that have been used for neutral dust agglomeration (e.g. hit and stick, ballistic motion[6]), wet granulate aggregation (sticking due to capillary bridges and ballistic motion[16]), colloidal aggregation (van der Waals and repulsions[34]), and hit and stick agglomeration in Brownian particles under frictional drag[35]. 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[25] tends to resemble a DLPCA-like behavior where the smaller size aggregates are larger in number, in contrast to a BCCA-like model where the size distribution is typically bell-shaped[6]. On the other hand, the monomer motion is found to be highly sub-diffusive[25] in agreement with the BCCA model. In addition, the Coulombic interactions will cause considerable deviations from the short-ranged or ballistic propagation typical of the BCCA or DLPCA models. We find that the long-range forces due to a bipolar charge distribution lead to the value of intermediate between the above two aggregation models, indicated by dashed lines in Fig. 5(d).

Interplay between fractals and mesoscopic flow

Apart from the long-range 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 modified 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 quasi-monodisperse assumption typically used in cluster-cluster aggregation models. The quasi-monodispersity however neglects the morphology and surface irregularities of the colliding aggregates. In practice, collisions between two aggregates with large size difference are also possible. Additionally, individual monomers might collide with large aggregates and result in fragmentation events. Our kinetic description neglects these specific scenarios. It is important to note that the main reason for the mismatch between kinetic theory and MD results for the density evolution is our neglect of size differences among aggregates in the kinetic theory.

Secondly, granular gases are characterized by the emergence of an advective flow[27,36] which we find, in the present case, induces the non-monotonicity in the temporal evolution of the number density [Fig. 3(c)]. Due to the advective flow, aggregates which are weakly connected are prone to fragmentation. This results in an intermediate regime where the concentration of aggregates increases instead of decreasing.

Excluding the two above mechanisms explains the slight deviation of our quasi-monodisperse Boltzmann theory from the non-monotonic behavior of after the crossover to aggregative collapse.

Discussion

We have derived the rate equations for the evolution of the number density , granular temperature , and charge variance of the cluster population in a collisionally charged, aggregating granular gas. In contrast to well-known Smoluchowski-type equations, we have explicitly coupled to the decay of and charge variance. We have compared the results with three-dimensional 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 self-similarity, 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 granulates with ion transfer mechanism[37,38], dissipative cell or active particle collections under long-range hydrodynamic and electrostatic effects[39,40], and charged ice-ice collisions[41].

Methods

Kinetics and modified collision integral

After obtaining the number of direct collisions  in Eq. (2) in the main text, let us consider the number of particles per unit spatial volume having initial velocity-charge values and in the intervals and which, post-collision, enter the intervals and in time Δt and thus the net increase of number of particles per unit time and volume is . We can relate the primed velocities to the unprimed viawhere is the Jacobian of the transformation for viscoelastic particles[27]. Here is a material constant. To obtain the transformation , we consider the ratio of relative charges after and before the collisionand in addition we impose charge conservation during collisions

The above two relations finally provide the transformationwhere, for example, for a constant r. This means that the differential charge-space volume element shrinks or expands by a factor of r/2. In general, the charge transfer may depend on pre-collision velocities and charges, and the expressions of r and can be quite complicated. The charge exchange depends on myriad factors such as size, composition, and crystalline properties. Incorporating the above phase-space volume transformations due to collisions, the net change ΔN of number of particles per unit phase-space volume and in time Δt readswhere we assume that the differential cross-section and the contact condition specified by are the same for direct and inverse collisions. Finally, dividing by Δt, and integrating over all incoming particle velocities and charges from all directions in the limit , we obtain a formal expression for the collision integral

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 for certain conditions on the relative velocity , and by considering the charge transferred to particle equal to the charge on particle [Eqs (25–30) below]. In Icoll, 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 long-range effect is incorporated via collision cross-section [Eq. (38)].

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.

It can be shown thatwhere and is the change of during the collision between pair i, j, and the prime denotes a post collision value Icoll. 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), respectivelywhere is the mean charge on the colliding pair. At this point we differentiate the restitutive or dissipative collisions from aggregative ones by splitting aswhere

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 non-zero , 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 (). Also represented schematically in Fig. 2 in the main text.

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.

For the granular temperature,where we take the limit for the aggregation. The change in the charge variance is obtained aswhere equals the charge transferred to particle i during its collision with particle . For the charge transfer, based on seminal experiments[10], we considerwhich is also obtainable if the charge transferred is considered proportional to the maximum contact area during the course of a collision[43]. In the present simulations, the value of the exponent after the experiments[10,42], where a power law dependence was found of the charge exchange on the relative collision energy, when silica particles were impacted on a surface. In the MD simulations, we have verified that while moderate changes of around the experimental value might produce changes in the charge buildup rate of the gas ( in Fig. 3(b)), the dynamics of charge in the aggregation phase are unaffected ( in Fig. 3(b)). Once the magnitude of collision velocities reduces with reducing granular temperature, the value of has negligible influence on the decay rate of , because it is primarily controlled by aggregation and merging of charges. The power law in Fig. 3(b) during the aggregation time regime remains unaffected.

Using Eq. (28) in Eq. (27), we  obtain

For aggregation, the charge tranferred to particle equals the charge on the merging particle , i.e., , which gives

Putting Eqs (25, 26, 29 and 30) in (23), and then (23) in (19), the resulting integrals are solved, assuming the statistical independence of charge-velocity distribution function, i.e., , and assuming that their scaled form remains Gaussian. In addition to the charge exchange, the coefficient of restitution is taken as velocity dependent[30,44], i.e.while the long-range effects due to Coulomb interactions are taken into account by the change in collision cross section [Eq. (38)]. After integration we obtain Eqs (6–8) in the main text. The functions g in Eqs (6–8) have the formswhere the coefficients , , , , are functions of [Table 1].

Derivation of the hydrodynamic equations (6–8)

To solve the integrals in Eq. (19) for different Δ from Eqs (25–30), we assume that the normalized velocity as well as charge distribution of the aggregating particle population at any time remain Gaussian, and the two are uncorrelated, i.e.

In Fig. 6, we show the charge distribution of monomers obtained from typical simulation runs, which is essentially the distribution before initiation of the aggregation. In Eq. (37) above, we assume that although granular temperature T and charge variance do change with time due to restitution and aggregation, the shape of the scaled distribution remains close to Gaussian, and the increase of size/decrease of number of particles due to aggregation process does not alter the shape of the scaled distribution.

Figure 6

The scaled charge distribution of individual particles (monomers) obtained from typical MD simulation runs (dots) in the aggregated granular gas. The solid line is a Gaussian fit. Here .

The attractive or repulsive long-range effects are emulated through an effective differential cross-section for a binary collision, which changes depending upon the sign and magnitude of charges on the particle pair i, j and their relative velocity, according towhere is the differential cross-section per unit solid angle . The expression is independent of n, and thus the total cross section is = = . For neutral particles, , and thus . For (repulsive encounters), , while for (attractive encounters), . Thus Eq. (38) is a linear adjustment to the neutral cross-section and is reasonable approximation in case of small angle scattering[45]. The possibility of negative cross-section for repulsive encounters is eliminated by the function  .

Below we explain the solution procedure for the restitutive, as well as aggregative, part of the equation for . Similar procedure can be followed for the equations for and .

Plugging Eq. (26) into Eq. (23) and then the resulting equation into Eq. (19), we findwhich, after using the above Eqs (31 and 38), and separating the integrals over n, v and q, reads aswhereandwhere in the aggregative part, we have set , and selects only the attractive encounters against low velocities selected by -- the charge-velocity combination which leads to aggregation.

Solution for the restitutive part

The solution for the parts , are as follows.

Using Eqs (37 and 46) from the above text, reads as

To perform the integration over the relative velocity , the following transformations are made: (i) , and (ii) , which also results in . Incorporating these transformations, we obtain

Notice the lower limit on relative velocities,is also altered due to the transformation . The integral giveswhile the integral is solved aswhere . Putting , in , and finally in the restitutive part of Eq. (40), and integrating over , , we obtain

Notice that if , we recover the classical Haff’s law for viscoelastic granular gas.

Solution for the aggregative part

The solution for the parts , are as follows.

Using Eqs (37 and 53) from the above text, and again using the variable transformations (i) , (ii) , (iii) , the integral in (42) reduces to

Here notice that now the relative velocity limits are from to , the condition for aggregation selected by . Finally putting (54) into aggregative part of (40) and integrating over q, q, we obtain

Notice that after integrating from to , the integration over q, q is to be broken into the sum of two parts, one over , , plus a second integral over , , to satisfy the aggregative condition set by .

A similar procedure is followed for , and using corresponding Δ (See Supplementary Mathematica script). Finally, the key constraint to be noted is that the solutions of the integrals in the rate equation for T are real valued for , while in equations for and , they are real valued for . The MD simulations confirm that these constraints put a physical limit during the aggregation phase.

Granular MD simulations

The equation of motion of the formis solved for each particle with a setup of periodic boundary conditions in a cubic box of size , where d0 is the monomer diameter and L is the non-dimensional system length [Tables 2 and 3]. Here b is a vector of integers representing the periodic replicas of the system in each Cartesian direction. The symbol ′ indicates that if to avoid Coulomb interaction of particles with themselves. The non-dimensional numbers in the above equation are , and , with α and β being viscoelatic material constants. For practicality, we select the reference length , time reference , velocity reference , and charge reference such that the elastic force strength , dissipative force strength /10, and Coulomb force strength is varied across [Tables 2 and 3]. The effect of dissipation compared to elastic forces is extensively studied for neutral systems[27]. We have reported on the variation of Coulomb strength compared to dissipation and elastic forces in[25]. The above order of magnitudes of , , and also help to attain an early clustering in non-dimensional time units in a finite size () neutral granular gas system [see[25] for more details]. Also , , n is the unit vector pointing from center of in-contact neighbor j towards the center of particle i, while r is the distance vector pointing from particle k towards the center of particle i.

The equation of charge on particle i may be written aswith I being the charge-exchange currents from colliding neighbors j during the course of collision. For any contact neighbor j, we approximate its integrated value over the time  interval by Eq. (28), i.e.

See[25] for more details. The long-range Coulomb forces for a setup with periodic boundary conditions in Eq. (56) 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 highly parallelized and optimized on graphics processing unit (GPU). In our simulations, the total computing time to reach non-dimensional simulation time ~103 for a typical simulation with monomers , including the long-range electrostatic forces, is of the order of weeks. See[25] for more details.

Comparison of individual T, 〈δq2〉, and n profiles, and emergence of advective flow using Mach number

In Fig. 7, we decompose the theoretical comparison of , presented in the main text, into individual comparisons of , , and profiles for a typical simulation run. The difference between the kinetic theory with and without aggregation is also emphasized.

Figure 7

Individual comparisons of (a) temperature of cluster population , (b) charge variance of cluster population , and (c) number density of cluster population , with the theoretical predictions. The granular temperature, charge variance and average size of the cluster population during aggregation evolve in such a manner that their non-dimensional combination (main text). The number density evolution, however, is highly dynamic and exhibits a non-monotonic behavior due to emergence of mesoscopic flow, marked by a transition in the rate of increase of Mach number Ma, shown in (c) inset. The numbers in (c) inset are the power law exponents of the time dependence of Ma. The short-hands res and agg denote restitution and aggregation respectively.

It is noticeable that the granular temperature of the aggregates closely follows Haff’s law, and is confirmed by theory, notwithstanding the presence of long-range effects and intricate aggregation and annihilation events. If only the restitutive terms of the hydrodynamic equations are considered (dashed line), the theory predicts that drops at a slower rate at long times. Furthermore, the charge variance in this case saturates. The number density in the absence of aggregation is, of course, invariant. If the aggregation dynamics is augmented, the simulation results are closely predicted by the theory.

In the MD simulations, the decay of during the aggregation phase closely agrees with the theory, even though we observe that the charge exchange in the simulations leads to a symmetric but non-Gaussian charge distribution among monomers during the initial restitution phase [Fig. 1].

The number density evolution, however, is highly dynamic and exhibits a non-monotonic behavior due to fragmentation events caused by the emergence of mesoscopic flow (see Suppl. Movie). The theory predicts the decay of cluster density only in an average sense [Fig. 7(c) inset]. To quantify the emergence of the mesoscopic flow, we calculate the Mach numberwhere V is the local mesoscopic velocity, and is the thermal velocity. To compute this, we divide the system volume into equal sized cubic boxes. The advective velocity in j box is computed as , where index i runs over all the monomers in that box j. We take the square of this velocity (we square the velocity of the box before summing it up over all the boxes, otherwise the advective velocities of the boxes might cancel each other even in the presence of a directed motion/advective flow, for example in a lattice of vortices) as , take its square root, and normalize it with the thermal velocity . Thus according to our definition, Ma is a global measure of the magnitude of directed/advective motion. The time evolution of Ma is shown in Fig. 7(c) (inset), which indicates that Ma grows at a higher rate when the number density evolution becomes non-monotonic (), indicating an intricate aggregation and fragmentation dynamics, and the generation of mesoscopic flow.

Reference scales and laboratory relevance of present results

Typically non-Brownian growth in planetary dust becomes dominant for monomer sizes near or above several μm[46] and the growth barrier problem[12] starts to arise near . The mass of silica particles in this range of sizes is . If the particles are initially agitated with velocities , the time scale reference to convert our simulation time to laboratory time is . Thus in our results the growth over 104 units of non-dimensional time approximately implies growth over ~10 s. The average size of aggregates in the growth period grows approximately by one order of magnitude (e.g. the growth of the largest cluster is from ≈2 mm to ≈7 cm in ≈10 s time for particles of such size and mass, and for initial monomer filling fraction of ) [Fig. 8].

Figure 8

Growth of the average cluster size, and the size of the largest cluster.

Statistical and analytical methods

The following method codes are used:

An in-house MATLAB program to solve Eqs (6–8).

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

Mathematica instructions to solve kinetic integrals

Movie of aggregating charged granular gas