Electrification in granular gases leads to constrained fractal growth

TL;DR

This paper develops a hydrodynamic model for charged granular gases, revealing how electrostatic forces constrain cluster growth and lead to fractal aggregation, supported by theory and molecular dynamics simulations.

Contribution

It introduces a novel theoretical framework for analyzing charge-driven aggregation in granular gases, highlighting a universal energy ratio constraint and fractal growth behavior.

Findings

Electrostatic forces constrain cluster energy ratios to be less than or equal to one.

Aggregation proceeds in a fractal manner influenced by Coulomb interactions.

The model is validated by molecular dynamics simulations.

Abstract

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 system obeys a physical constraint of nearly constant dimensionless ratio of characteristic electrostatic to kinetic energy $\mathcal{B}(t)\le 1$. This constraint on the collective evolution of charged clusters is confirmed both by the theory and the 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 mechanism, 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.

Figures

Tables

Table 1. Table 1: Expressions of the coefficients in Eq. (5)-(7) in the main text. Here m 𝑚 m and d 𝑑 d are the mass and size of the aggregates. The material constant C ϵ subscript 𝐶 italic-ϵ C_{\epsilon} is from Eq. ( 30 ) and influence the viscoelastic properties of the particles, while C Δ ​ q subscript 𝐶 Δ 𝑞 C_{\Delta q} and η 𝜂 \eta are from Eq. ( 27 ) and influence the charge buildup. Other notations are as described in the main text.

Coefficient	Expression
$ℬ$	$k_e⟨\delta q^2⟩/\left(Td\right)$
$l$	$\sqrt{4-{ℬ}^2}$
$l_1$	$\sqrt{1-{ℬ}^2}$
$C_{res}^T$	${\mathrm{80\hspace{0.33em}2}}^{1/5}{d}^2{C}_{ϵ}m/\left(21\sqrt{\pi }{m}^{8/5}\right)$
$C_{agg}^T$	$d^2/\left(8\sqrt{\pi m}\right)$
$C_{res}^q$	$1.4080d^2{C}_{\Delta q}^2\left(2+\eta \right){\left(4/m\right)}^{\eta +1/2}/\pi$
$C_{agg}^q$	$4d^2/\sqrt{\pi m}$
$C_{agg}^n$	$d^2/\left(2\sqrt{\pi m}\right)$
$a_1^T$	$ℬ\left(8-5{ℬ}^2\right)l+\pi \left(16-10{ℬ}^2+3{ℬ}^4\right)$
$a_2^T$	$-32+20{ℬ}^2-6{ℬ}^4$
$a_3^T$	$l\left(2\pi {\left({ℬ}^2-4\right)}^2+ℬ\left(-32+28{ℬ}^2+{ℬ}^4\right)\right)$
$a_4^T$	$-4\left(32-20{ℬ}^2+9{ℬ}^4\right)$
$a_5^T$	$lℬ\left(-32+28{ℬ}^2+{ℬ}^4\right)$
	$+2\pi \left({ℬ}^2\left(20-8l\right)+{ℬ}^4\left(-9+l\right)+16\left(-2+l\right)\right)$
$a_6^T$	$4\left(32-20{ℬ}^2+9{ℬ}^4\right)$
$a_1^q$	$2ℬl\left({ℬ}^2-1\right)-\pi \left(4-2{ℬ}^2+{ℬ}^4\right)$
$a_2^q$	$2\left(4-2{ℬ}^2+{ℬ}^4\right)$
$a_3^q$	$16l\sqrt{\pi }{l}_1{\left(4-5{ℬ}^2+{ℬ}^4\right)}^2$
$a_4^q$	$-\sqrt{\pi }ℬl(\pi ({ℬ}^2-4)(1+2{ℬ}^2)$
	$+2ℬl_1(-8-54{ℬ}^2+33{ℬ}^4+2{ℬ}^6))$
$a_5^q$	$2\sqrt{\pi }ℬl^5\left(1+2{ℬ}^2\right)$
$a_6^q$	$16\sqrt{\pi }ℬl_1^5\left(4+5{ℬ}^2\right)$
$a_7^q$	$8l\sqrt{\pi }{l}_1{\left(4-5{ℬ}^2+{ℬ}^4\right)}^2$
$a_8^q$	$-\sqrt{\pi }ℬl^5\left(1+2{ℬ}^2\right)$
$a_9^q$	$\sqrt{\pi }ℬl_1ℬl\left(8+54{ℬ}^2-33{ℬ}^4-2{ℬ}^6\right)$
$a_{10}^q$	$8\sqrt{\pi }ℬl_1{\left({ℬ}^2-1\right)}^2\left(4+5{ℬ}^2\right)$
$a_1^n$	$-{ℬ}^2{l}^3$
$a_2^n$	$l_{1}lℬ\left(-4+7{ℬ}^2\right)+l_1\pi \left({ℬ}^2-1\right)\left(8-4l+{ℬ}^2\left(-6+l\right)\right)$
$a_3^n$	$4l_1\left(4-7{ℬ}^2+3{ℬ}^4\right)$
$a_4^n$	$-{ℬ}^2{l}^3$
$a_5^n$	$l_{1}l\left(ℬ\left(-4+7{ℬ}^2\right)+\pi \left(4-5{ℬ}^2+{ℬ}^4\right)\right)$
$a_6^n$	$-4l_1\left(4-7{ℬ}^2+3{ℬ}^4\right)$

Table 2. Table 2: Simulation parameters.

Parameter	Expression	Value
$N_{\text{mon}}$	No. of monomers	$50016$
$L$	System length	$70$
$\varphi$	Monomer filling fraction	0.013, 0.030, 0.076, 0.137
$ℰ$	$\frac{\alpha l_{\text{ref}}^{3/2}{t}_{\text{ref}}}{m_{\text{ref}}{v}_{\text{ref}}}$	278
$𝒟$	$\frac{\beta l_{\text{ref}}^{1/2}{t}_{\text{ref}}}{m_{\text{ref}}}$	27.8
$𝒦$	$\frac{k_e{q}_{\text{ref}}^2{t}_{\text{ref}}}{m_{\text{ref}}{v}_{\text{ref}}{l}_{\text{ref}}^2}$	0.4, 1.0, 5.0

Table 3. Table 3: Reference values in SI units for non-dimensionalization.

Reference scale	Description	Value
$m_{\text{ref}}$	Particle mass reference	$1.52\times {10}^{-4}$ $\mathrm{kg}$
$l_{\text{ref}}=d_0$	Length reference=monomer diameter	$4.78\times {10}^{-3}$ $m$
$v_{\text{ref}}$	Velocity reference	$1.0$ ${\mathrm{ms}}^{-1}$
$t_{\text{ref}}=d_0/v_{\text{ref}}$	Time reference reference	$4.78\times {10}^{-3}$ $s$
$T_{\text{ref}}=\frac{1}{3}{m}_{\text{ref}}{v}_{\text{ref}}^2$	Temperature reference	$0.51\times {10}^{-4}$ $\text{kg}{\text{m}}^2{\text{s}}^{-2}$
$k_e$	Coulomb’s constant	$8.98\times {10}^9$ $\text{N}{\text{m}}^2{\text{C}}^{-2}$
$\alpha$	Elastic constant of particles	$2.67\times {10}^4$ $\text{kg}{\text{m}}^{-1/2}{\text{s}}^{-2}$
$\beta$	Viscous constant of particles	$1.28\times {10}^1$ $\text{kg}{\text{m}}^{-1/2}{\text{s}}^{-1}$
$q_{\text{ref}}$	Charge reference	$2.0\times {10}^{-8}$ C - $5.6\times {10}^{-9}$ C

Equations

$$\frac{\partial f}{\partial t}=I_{\mathrm{coll}},$$

$$N_c^{-}=f_{i}d{\mathbf{v}}_{i}d{q}_i{f}_{j}d{\mathbf{v}}_{j}d{q}_j|{\mathbf{v}}_{ij}\cdot \mathbf{n}|\Theta (-{\mathbf{v}}_{ij}\cdot \mathbf{n}){\Theta }_{q}d\sigma \Delta t,$$

$$N_c^{+}=f_i^{\prime \prime }d{\mathbf{v}}_i^{\prime \prime }d{q}_i^{\prime \prime }{f}_j^{\prime \prime }d{\mathbf{v}}_j^{\prime \prime }d{q}_j^{\prime \prime }|{\mathbf{v}}_{ij}^{\prime \prime }\cdot \mathbf{n}|\Theta (-{\mathbf{v}}_{ij}^{\prime \prime }\cdot \mathbf{n}){\Theta }_q^{\prime \prime }d{\sigma }^{\prime \prime }\Delta t.$$

$$\Delta N_c=$$

$$\times \Theta (-{\mathbf{v}}_{ij}\cdot \mathbf{n})d{\mathbf{v}}_{j}d{q}_{j}d\sigma {\Theta }_q\Delta t,$$

$$I_{\mathrm{coll}}=$$

$$\times \Theta (-{\mathbf{v}}_{ij}\cdot \mathbf{n})d{\mathbf{v}}_{j}d{q}_{j}d\sigma {\Theta }_q.$$

$$\frac{\partial n}{\partial t}$$

$$\frac{3}{2}\frac{\partial T}{\partial t}=-n^2{T}^{\frac{8}{5}}{g}_2(\mathcal{B},C_{res}^T)+n^2{T}^{\frac{3}{2}}{g}_3(\mathcal{B},C_{agg}^T),$$

$$\frac{\partial ⟨\delta q^2⟩}{\partial t}=n^2{T}^{(\eta +\frac{1}{2})}{g}_4(\mathcal{B},C_{res}^q)-n^2⟨\delta q^2⟩T^{\frac{1}{2}}{g}_5(\mathcal{B},C_{agg}^q),$$

$$\mathcal{B}(t)\equiv k_e\frac{⟨\delta q^2⟩(t)}{T(t)d(t)},$$

$$d(t)={\left[\frac{6M}{\pi n(t)V{\rho }_p}\right]}^{\frac{1}{3}},$$

$$\mathcal{B}(t)\le 1.$$

$$N_c^{+}=f_i^{\prime \prime }d{\mathbf{v}}_i^{\prime \prime }d{q}_i^{\prime \prime }{f}_j^{\prime \prime }d{\mathbf{v}}_j^{\prime \prime }d{q}_j^{\prime \prime }|{\mathbf{v}}_{ij}^{\prime \prime }\cdot \mathbf{n}|\Theta (-{\mathbf{v}}_{ij}^{\prime \prime }\cdot \mathbf{n}){\Theta }_q^{\prime \prime }d{\sigma }^{\prime \prime }\Delta t,$$

$$d{\mathbf{v}}_i^{\prime \prime }d{\mathbf{v}}_j^{\prime \prime }={\mathbf{J}}_{ij}^{v_{ij}}d{\mathbf{v}}_{i}d{\mathbf{v}}_j,$$

$$r=\frac{q_i-q_j}{q_i^{\prime \prime }-q_j^{\prime \prime }},$$

$$q_i+q_j=q_i^{\prime \prime }+q_j^{\prime \prime }.$$

$$d{q}_i^{\prime \prime }d{q}_j^{\prime \prime }={\mathbf{J}}_{ij}^{q}d{q}_{i}d{q}_j,$$

$$\Delta N_c=$$

$$I_{\mathrm{coll}}=$$

$$\frac{\partial ⟨\psi ⟩}{\partial t}=\int d{\mathbf{v}}_{i}d{q}_i{\psi }_i\frac{\partial f_i}{\partial t}=\int d{\mathbf{v}}_{i}d{q}_i{\psi }_i{I}_{\mathrm{coll}}.$$

$$\frac{\partial ⟨\psi ⟩}{\partial t}$$

$$=\frac{1}{2}\int d{\mathbf{v}}_{i}d{\mathbf{v}}_{j}d{q}_{i}d{q}_{j}d\sigma f_i{f}_j|{\mathbf{v}}_{ij}\cdot \mathbf{n}|\Theta (-{\mathbf{v}}_{ij}\cdot \mathbf{n}){\Theta }_q\Delta [{\psi }_i+{\psi }_j]$$

$$=\int d{\mathbf{v}}_{i}d{\mathbf{v}}_{j}d{q}_{i}d{q}_{j}d\sigma f_i{f}_j|{\mathbf{v}}_{ij}\cdot \mathbf{n}|\Theta (-{\mathbf{v}}_{ij}\cdot \mathbf{n}){\Theta }_q\Delta [{\psi }_i],$$

$$(i)\psi$$

$$(ii)\psi$$

$$(iii)\psi$$

$${\Theta }_q\Delta [{\psi }_i]$$

$$+{\Delta }^{agg}[{\psi }_i]\Theta (-q_i{q}_j)\Theta \left(\sqrt{b}-v_{ij}\right),$$

$$\sqrt{b}\equiv \sqrt{\frac{2k_e|q_i{q}_j|}{md}}.$$

$${\Delta }_n^{res}[{\psi }_i+{\psi }_j]$$

$${\Delta }_n^{agg}[{\psi }_i+{\psi }_j{]}_n$$

$${\Delta }_T^{res}[{\psi }_i+{\psi }_j]$$

$${\Delta }_T^{agg}[{\psi }_i+{\psi }_j]$$

$${\Delta }_q^{res}[{\psi }_i]$$

Full text

I 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 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, 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 Jungmann et al. (2018); Lee et al. (2015); Yoshimatsu et al. (2017); Poppe et al. (2000); Haeberle et al. (2018). However, technical difficulties plague even careful experiments and often impede their unambiguous interpretation Spahn and Seiß (2015). 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 Spahn and Seiß (2015). Despite considerable effort Ivlev et al. (2002); Dammer and Wolf (2004); Müller and Luding (2011); Ulrich et al. (2009); Brilliantov et al. (2018); Liu and Hrenya (2018); Takada et al. (2017); Kolehmainen et al. (2018) 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 velocities $\epsilon(v)$, 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 $\langle\delta q^{2}\rangle$ of the aggregates under the assumptions of homogeneous and quasi-monodisperse aggregation. 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 characteristic electrostatic to thermal energy.

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 macroscopic flow.

II 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 Singh and Mazza (2018). 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 supplementary for MD). In the following we establish a modified Boltzmann approach for this intricate dynamics of aggregation process, which predicts novel physical limits.

We consider the single particle probability distribution function $f=f(\mathbf{r},t;\mathbf{v},q,d)$, where the particle velocity $\mathbf{v}$, charge $q$, position $\mathbf{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. $f=f(t;\mathbf{v},q)$ only, and its time evolution is given by the simplified Boltzmann equation Pitaevskii and Lifshitz (2012); Brilliantov and Pöschel (2010)

[TABLE]

valid at any time instant $t$. Here we define $I_{\mathrm{coll}}$ 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.

Let us consider contact collisions of particles $i$ and $j$ with pre-collision velocity-charge values $(\mathbf{v}_{i},q_{i})$ and $(\mathbf{v}_{j},q_{j})$, respectively. In the ensemble picture, particle collisions will change $f(t;\mathbf{v},q)$ in the infinitesimal phase-space volumes $d\mathbf{v}_{i}dq_{i}$ and $d\mathbf{v}_{j}dq_{j}$, centered around $(\mathbf{v}_{i},q_{i})$ and $(\mathbf{v}_{j},q_{j})$, respectively. The number of direct collisions $N_{c}^{-}$ per unit spatial volume which lead to loss of particles from the intervals $d\mathbf{v}_{i}dq_{i}$ and $d\mathbf{v}_{j}dq_{j}$ in time $\Delta t$ are

[TABLE]

where $\mathbf{v}_{ij}\equiv\mathbf{v}_{i}-\mathbf{v}_{j}$, $\mathbf{n}$ is the unit vector at collision pointing from the center of particle $i$ towards particle $j$, and $d\sigma$ is the differential collision cross-section. The Heaviside step function $\Theta(-\mathbf{v}_{ij}\cdot\mathbf{n})$ selects particles coming towards $i$, while we use $\Theta_{q}\equiv\Theta\left(\frac{1}{2}mv_{ij}^{2}-\frac{k_{e}q_{i}q_{j}}{d}\right)$ to ensure that a contact with an approaching particle takes place only when the Coulomb energy barrier can be overcome, where $k_{e}=1/(4\pi\epsilon_{0})$, $\epsilon_{0}=8.854\times 10^{-12}$ F m*-1* is the vacuum permittivity, and $d$ is the particle diameter at time $t$. If the interaction is repulsive, ${k_{e}q_{i}q_{j}}/{d}$ is positive, and $\Theta_{q}=1$ only if $\frac{1}{2}mv_{ij}^{2}>{k_{e}q_{i}q_{j}}/{d}$. In case of attractive interaction, ${k_{e}q_{i}q_{j}}/{d}$ is negative and thus $\Theta_{q}=1$ always. Essentially, $\Theta_{q}$ filters repulsive interactions which do not lead to a physical contact between particles.

Consider now particles with initial velocity-charge values $(\mathbf{v}_{i}^{\prime\prime},q_{i}^{\prime\prime})$ and $(\mathbf{v}_{j}^{\prime\prime},q_{j}^{\prime\prime})$ in the intervals $d\mathbf{v}_{i}^{\prime\prime}dq_{i}^{\prime\prime},\;d\mathbf{v}_{j}^{\prime\prime}dq_{j}^{\prime\prime}$. The number of particles $N_{c}^{+}$ per unit volume which, post-collision, enter the interval $d\mathbf{v}_{i}dq_{i}$ and $d\mathbf{v}_{j}dq_{j}$ in time $\Delta t$ is

[TABLE]

The net change $\Delta N_{c}\equiv N_{c}^{+}-N_{c}^{-}$ of number of particles in time $\Delta t$ per unit volume, then reads (see Supplementary)

[TABLE]

where $\mathbf{J}^{v_{ij}}_{ij}$ and $\mathbf{J}^{q}_{ij}$ are the Jacobians of the transformations for $d\mathbf{v}_{i}^{\prime\prime}d\mathbf{v}_{j}^{\prime\prime}\rightarrow d\mathbf{v}_{i}d\mathbf{v}_{j}$ and $dq_{i}^{\prime\prime}dq_{j}^{\prime\prime}\rightarrow dq_{i}dq_{j}$, 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 $\Delta t$ and taking the limit $\Delta t\to 0$, we obtain the formal expression for the modified collision integral

[TABLE]

Here we assume that the differential collision cross section and the contact condition specified by $\Theta_{q}$ retain their form for direct and inverse collisions. The particle encounters which do not lead to a physical contact have been excluded using $\Theta_{q}$. While taking moments of $I_{\mathrm{coll}}$, a fraction of those contact collisions that lead to aggregation is accounted for by taking the limit $\epsilon=0$ for certain conditions on the relative velocity $v_{ij}$, and by considering the charge transferred to particle $i$ equal to the charge on particle $j$ (see Supplementary). In $I_{\mathrm{coll}}$, 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.

After setting up the collision integral, we derive the macroscopic changes of number density $n$, granular temperature $T$, and the charge variance $\langle\delta q^{2}\rangle$ by taking the moments of the Boltzmann equation (see Supplementary). 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 $\langle\delta q\rangle$ is zero. The next choice is thus $\langle\delta q^{2}\rangle$. 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 (see Supplementary). After integration we find the governing equations

[TABLE]

[TABLE]

[TABLE]

which are coupled via a time-dependent dimensionless ratio

[TABLE]

between charge variance, granular temperature, and aggregate size. The terms $g_{k}$ are time-dependent functions of $\mathcal{B}$ and material constants $C_{res},\;C_{agg}$ (Supplementary Table I). We term the ratio $\mathcal{B}$ as Bjerrum number. In Eq. (8), $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 $\rho_{p}$ to be constant, which fixes the relation between particle size $d$ and particle number density $n$, according to

[TABLE]

and closes the equation set (5)-(7). The above set of equations is consistent with the classical Haff’s law in the absence of collisional charging. In this limit $\langle\delta q^{2}\rangle=0$, $\mathcal{B}=0$, and we obtain $\frac{\partial n}{\partial t}=0$, $\frac{\partial\langle\delta q^{2}\rangle}{\partial t}=0$ and $\frac{3}{2}\frac{\partial T}{\partial t}=-T^{8/5}\left[\frac{\pi C^{T}_{res}}{2}\right]$, whose solution is Haff’s law.

III Results and Discussion

The time evolution of $T$, $\langle\delta q^{2}\rangle$ and $n$ of aggregate population from MD simulations is shown in Fig. 3. The aggregate temperature from simulations is extracted as $\frac{3}{2}T=\frac{1}{N_{\mathrm{agg}}}\sum_{i}\frac{1}{2}m_{i}[(\mathbf{v}_{i}-\mathbf{V})^{2}]$. Notice that $\mathbf{v}_{i}$ and $m_{i}$ are center of mass velocity, and mass of the $i^{th}$ aggregate respectively, and should not be confused with monomer velocities and masses. $\mathbf{V}$ is the local advective velocity in the neighborhood of $i^{th}$ aggregate. Similarly, $\langle\delta q^{2}\rangle=\frac{1}{N_{\mathrm{agg}}}\sum_{i}(q_{i}-\langle q\rangle)^{2}$, and $n=N_{\mathrm{agg}}/V$, where $N_{\mathrm{agg}}$ is the total number of aggregates and $V$ is the system volume.

Initially ($tv_{\textrm{ref}}/l_{\textrm{ref}}<10^{2}$), the relative collision velocities $v_{ij}$ remain larger than the time varying threshold $\sqrt{b}\equiv\sqrt{\frac{k_{e}|q_{i}q_{j}|}{2Td}}$ (see Supplementary for the treatment of threshold $b$). 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 $\langle\delta q^{2}\rangle$ [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$, $T$ and $\langle\delta q^{2}\rangle$ can be collected into evolution of $\mathcal{B}$, shown in Fig. 3(d). The Bjerrum number $\mathcal{B}$ 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_{ij}$ approach the threshold $\sqrt{b}$, $\mathcal{B}\to 1$. 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 next section and in the Supplementary. These results are robust under variation of the initial monomer filling fraction $\phi$ and the charge strength $\mathcal{K}$ [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 the combination $\mathcal{B}$, and its comparison with the solution of Eq. (5)-(7) is highlighted. We solve the full kinetic equations Eq. (5)-(7) including the aggregation kinetics [Fig. 4 (solid line)]. We also solve Eq. (5)-(7) 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, $\mathcal{B}$ 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

[TABLE]

The upper physical limit predicted in the theory, $\mathcal{B}(t)<1$, is endorsed by the granular MD simulations under moderate variation of $\phi$ and $\mathcal{K}$. It is notable that at later times, the limit $\mathcal{B}\leq 1$ allows the right hand sides of the equations for number density, temperature and charge variance [Eq. (5)-(7)] 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 thermal motion of monomers, if any, inside the clusters is significantly decreased, and thus collisional charge-exchange is expected to be negligible. Also, in our kinetic theory description, the charge transferred to particle $i$ during aggregation with particle $j$ is equal to the charge on particle $j$, as they merge into one single particle [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 $\mathcal{B}(t)$. The constraint $\mathcal{B}\leq 1$ is robust in the theory, while the granular MD simulations suggest $\mathcal{B}(t)<1$ and confirm the upper limit of $\mathcal{B}(t)$. It is also intriguing that the temperature of the cluster population still closely follows Haff’s law for different $\phi$ and $\mathcal{K}$, 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. It initially begins to decrease during small aggregate formations due to low velocity attractive monomer encounters. In an intermediate time regime, the emergence of macroscopic particle fluxes trigger fragmentation events and the aggregate numbers increase. We quantify the emergence of macroscopic flow using the Mach number (see Supplementary 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 $\mathcal{B}(t)$ after the aggregative collapse ($tv_{\mathrm{ref}}/l_{\mathrm{ref}}>10^{2}$) [Fig. 3 and 4]. The distribution of clusters neglected in the homogeneous and quasi-monodisperse aggregation picture adds to the complexity of $\mathcal{B}$’s evolution. However, the theory clearly predicts the growth of $\mathcal{B}$ and selects an upper limit. 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.

III.1 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 $\langle D_{f}\rangle$ Mandelbrot (1977); Jullien (1987) of cluster population from the scaling relation $m\sim R_{g}^{\langle D_{f}\rangle}$ between cluster masses $m$, and radii of gyration $R_{g}=[\frac{1}{N_{\textrm{mon}}}\sum_{i}(\mathbf{r}_{i}-\mathbf{r}_{\textrm{mean}})^{2}]^{1/2}$, where the index $i$ runs over total number of monomers $N_{\textrm{mon}}$ in a given aggregate. Figure 5(a-c) show scatter plots for $R_{g}$ versus $m$ at different times, and Fig. 5(d) the time evolution of the exponent $\langle D_{f}\rangle$ for varying filling fraction. The average fractal dimension lies between the average values reported for the ballistic cluster-cluster aggregation (BCCA, $\langle D_{f}\rangle\simeq 1.94$) and the diffusion-limited particle-cluster aggregation (DLPCA, $\langle D_{f}\rangle\simeq 2.46$) Blum (2006); Smirnov (1990) models. In time, $\langle D_{f}\rangle$ 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 are popular models for aggregation that have been used for 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 $\langle D_{f}\rangle$ lies between the reported average values of $\langle D_{f}\rangle$ 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 Singh and Mazza (2018) 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 Blum (2006). On the other hand, the monomer motion is found to be highly sub-diffusive Singh and Mazza (2018) 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 $\langle D_{f}\rangle$ intermediate between the above two aggregation models, indicated by dashed lines in Fig. 5(d).

III.2 Interplay between fractals and macroscopic 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. Collisions between aggregates with large size difference, between aggregates and individual monomers, and the annihilation events are also simplified.

Secondly, granular gases are characterized by the emergence of a convective flow Hummel et al. (2016); Brilliantov and Pöschel (2010) which we find (see Supplementary and Media therein), in the present case, induces the non-monotonicity in the temporal evolution of the number density [Fig. 3(c)]. Due to the macroscopic 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 $\mathcal{B}(t)$ after the crossover to aggregative collapse.

IV Conclusions

We have derived the rate equations for the evolution of the number density $n$, granular temperature $T$, and charge variance $\langle\delta q^{2}\rangle$ of the cluster population in a charged, aggregating granular gas. In contrast to well-known Smoluchowski-type equations, we have explicitly coupled $n$ to the decay of $T$ 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 $\mathcal{B}(t)=k_{e}\langle\delta q^{2}\rangle/(Td)\leq 1$, 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 $T$ and $\langle\delta q^{2}\rangle$ 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 Lee et al. (2018); Zhang et al. (2015), dissipative cell or active particle collections under long-range hydrodynamic and electrostatic effects Yan et al. (2016); Friedl and Gilmour (2009), and charged ice-ice collisions Dash et al. (2001).

V Appendix

Kinetics and modified collision integral

After obtaining number of direct collisions in Eq. (2) in the main text, let us consider the number of particles $N_{c}^{+}$ per unit spatial volume having initial velocity-charge values $(\mathbf{v}_{i}^{\prime\prime},q_{i}^{\prime\prime})$ and $(\mathbf{v}_{j}^{\prime\prime},q_{j}^{\prime\prime})$ in the intervals $d\mathbf{v}_{i}^{\prime\prime}dq_{i}^{\prime\prime}$ and $d\mathbf{v}_{j}^{\prime\prime}dq_{j}^{\prime\prime}$ which, post-collision, enter the in the intervals $d\mathbf{v}_{i}dq_{i}$ and $d\mathbf{v}_{j}dq_{j}$ in time $\Delta t$ are

[TABLE]

and thus the net increase of number of particles per unit time and volume is $N_{c}^{+}-N_{c}^{-}$. We can relate the primed velocities to the unprimed via

[TABLE]

where $\mathbf{J}^{v_{ij}}_{ij}=1+(6/5)C_{\epsilon}v_{ij}^{1/5}+...$ is the Jacobian of the transformation for viscoelastic particles Brilliantov and Pöschel (2010). Here $C_{\epsilon}$ is a material constant. To obtain the transformation $d{q}_{i}^{\prime\prime}d{q}_{j}^{\prime\prime}\rightarrow d{q}_{i}d{q}_{j}$, we consider the ratio of relative charges after and before the collision

[TABLE]

and in addition we impose charge conservation during collisions

[TABLE]

The above two relations finally provide the transformation

[TABLE]

where, for example, $\mathbf{J}^{q}_{ij}=\frac{2}{r}$ for a constant $r$. This means that the differential charge-space volume element shrinks or expands by a factor of $r/2$. In general, for velocity and particle pre-charge dependent charge transfer, the expressions of $r$ and $\mathbf{J}^{q}_{ij}$ 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 phase-space volume transformations due to collisions, the net change $\Delta N_{c}$ of number of particles per unit phase-space volume and in time $\Delta t$ reads

[TABLE]

where we assume that the differential cross-section and the contact condition specified by $\Theta_{q}$ are the same for direct and inverse collisions. Finally, dividing by $\Delta t$, and integrating over all incoming particle velocities and charges from all directions in the limit $\Delta t\to 0$, we obtain a formal expression for the collision integral

[TABLE]

At this point the particle encounters which do not lead to a physical contact have been excluded using $\Theta_{q}$, however, collisions that lead to aggregation have not been explicitly accounted. We do so by taking the limit $\epsilon=0$ for certain conditions on the relative velocity $v_{ij}$, and by considering the charge transferred to particle $i$ equal to the charge on particle $j$ [Eq. (24)-(29) below]. In $I_{\mathrm{coll}}$, 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.

VI Splitting restitution and aggregation

The time rate of change of the average of a microscopic quantity $\psi(\mathbf{v}_{i},q_{i})$ is obtained by multiplying the Boltzmann equation for $f_{i}$ by $\psi_{i}$ and integrating over $\mathbf{v}_{i},q_{i}$, i.e.

[TABLE]

It can be shown that

[TABLE]

where $\Delta[\psi_{i}+\psi_{j}]=(\psi_{i}^{\prime}+\psi_{j}^{\prime}-\psi_{i}-\psi_{j})$ and $\Delta[\psi_{i}]=(\psi_{i}^{\prime}-\psi_{i})$ is the change of $\psi$ during the collision between pair $i,j$, and the prime denotes a post collision value. We note that the transformations in Eq. (12) and (15) are reversed back while integrating $I_{\mathrm{coll}}$ weighted with quantity of interest $\psi$. 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 $\langle\delta q\rangle$ is zero), respectively

[TABLE]

At this point we differentiate the restitutive or dissipative collisions from aggregative ones by splitting $\Theta_{q}\Delta[\psi_{i}]$ as

[TABLE]

where

[TABLE]

If $v_{ij}>\sqrt{b}$, the particles collide and separate after the collision irrespective of the sign of $q_{i}q_{j}$ (attractive or repulsive). This leads to dissipation of energy with finite non-zero $\epsilon=\epsilon(v_{ij})$, and charge exchange according to a specified rule. The aggregative part is zero in this case. If $v_{ij}<\sqrt{b}$ and $q_{i}q_{j}<0$ (attractive encounters at low velocities), the particles collide and aggregate with $\epsilon=0$, and with charge exchange to particle $i$ equal to $q_{j}$. If $v_{ij}<\sqrt{b}$ and $q_{i}q_{j}>0$ (repulsive encounters at low velocities), no physical contact takes place between the particles which leads to neither dissipation nor aggregation ($\Theta_{q}\Delta[\psi_{i}]=0$). Also represented schematically in Fig. (2) in the main text.

The expressions for $\Delta^{res}[\psi_{i}]$ and $\Delta^{agg}[\psi_{i}]$ are obtained as follows. The particle number does not change during a dissipative collision but reduces by one in an aggregative collision, i.e.

[TABLE]

For the granular temperature,

[TABLE]

where we take the limit $\epsilon=0$ for the aggregation. The change in the charge variance is obtained as

[TABLE]

where $(q_{i}^{\prime}-q_{i})$ equals the charge transferred to particle $i$ during its collision with particle $j$, and $q_{0}=\frac{q_{i}+q_{j}}{2}$ is the mean charge on the pair. For the charge transfer, based on seminal experiments Blum (2006), we consider

[TABLE]

which is also obtainable if charge transferred is considered proportional to the contact area during the course of collisions. Using Eq. (27) in Eq. (26), we get

[TABLE]

For aggregation, the charge tranferred to particle $i$ equals the charge on the merging particle $j$, i.e, $q_{i}^{\prime}-q_{i}=q_{j}$, which gives

[TABLE]

Putting Eq. (24),(25),(28), and (29) in (22), and then (22) in (18), the resulting integrals are solved, assuming the statistical independence of charge-velocity distribution function, i.e., $f(\mathbf{v},q)=f(\mathbf{v})f(q)$, and assuming that their scaled form remains Gaussian. In addition to the charge exchange, the coefficient of restitution is taken as velocity dependent, i.e

[TABLE]

while the long-range effects due to Coulomb interactions are taken into account by the change in collision cross section. After integration we obtain Eq. (5)-(7) in the main text. The functions $g_{k}$ in Eq. (5)-(7) have the forms

[TABLE]

where the coefficients $l$, $l_{1}$, $a^{T}_{k}$, $a^{q}_{k}$, $a^{n}_{k}$ are functions of $\mathcal{B}(t)$ [Table 1].

VII Derivation of the hydrodynamic equations (5)-(7)

To solve the integrals in Eq. (18) for different $\Delta[\psi_{i}]$ from Eq. (24)-(29), 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

[TABLE]

In Fig. 6, we show charge distribution of monomers obtained from typical simulation runs, which is essentially the distribution before initiation of the aggregation. In Eq. (36) above, we assume that although granular temperature $T$ and charge variance $\langle\delta q^{2}\rangle$ 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 scaled distribution shape.

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 to

[TABLE]

where $\left(\frac{d\sigma}{d\Omega}\right)$ is the differential cross-section per unit solid angle $d\Omega\equiv d\mathbf{n}$. The expression $\frac{d^{2}}{4}\left(1-\frac{2k_{e}q_{i}q_{j}}{dmv_{ij}^{2}}\right)$ is independent of $\mathbf{n}$, and thus the total cross section is $\sigma=\frac{d^{2}}{4}\left(1-\frac{2k_{e}q_{i}q_{j}}{dmv_{ij}^{2}}\right)\int_{0}^{\pi}d\phi d\theta\sin\theta=\pi{d^{2}}\left(1-\frac{2k_{e}q_{i}q_{j}}{dmv_{ij}^{2}}\right)$. For neutral particles, $q_{i}=q_{j}=0$, and thus $\sigma=\pi d^{2}$. For $q_{i}q_{j}>0$ (repulsive encounters), $\sigma<\pi d^{2}$, while for $q_{i}q_{j}<0$ (attractive encounters), $\sigma>\pi d^{2}$.

Below we explain the solution procedure for the restitutive, as well as aggregative, part of the equation for $\frac{\partial T}{\partial t}$. Similar procedure can then followed for the equations for $\frac{\partial n}{\partial t}$ and $\frac{\partial\langle\delta q^{2}\rangle}{\partial t}$.

Plugging Eq. (25) into Eq. (22) and then the resulting equation to Eq. (18), we find

[TABLE]

which, after using the above Eq. (30) and (37), and separating the integrals over $\mathbf{n},\mathbf{v}$ and $q$, reads as

[TABLE]

where

[TABLE]

and

[TABLE]

where in the aggregative part, we have set $\epsilon=0$, and $\Theta(-q_{i}q_{j})$ selects only the attractive encounters against low velocities selected by $\Theta(\sqrt{b}-v_{ij})$, the charge-velocity combination which leads to aggregation. Here

[TABLE]

VIII Solution for the restitutive part $\left(\frac{3}{2}\frac{\partial T}{\partial t}\right)_{res}$

The solution for the parts $I_{\mathbf{n}}^{res},I_{\mathbf{v}}^{res}$ are as follows.

[TABLE]

Using Eq. (45) and (36) from the above text, $I_{\mathbf{v}}^{res}$ reads as

[TABLE]

To perform the integration over the relative velocity $\mathbf{v}_{ij}$, the following transformations are made: (i) $\mathbf{w}_{ij}=\frac{2T}{m}(\mathbf{v}_{i}-\mathbf{v}_{j})$, and (ii) $\mathbf{w}_{c}=\frac{2T}{m}(\mathbf{v}_{i}+\mathbf{v}_{j})$, which also results in $d\mathbf{v}_{i}d\mathbf{v}_{j}=-\frac{1}{8}\left(\frac{m}{4T}\right)^{-3}d\mathbf{w}_{ij}d\mathbf{w}_{c}$. Incorporating these transformations, we obtain

[TABLE]

Notice the lower limit on relative velocities,

[TABLE]

is also altered due to the transformation $\mathbf{v}_{ij}\rightarrow\mathbf{w}_{ij}$. The integral $I_{w_{c}}$ gives

[TABLE]

while the integral $I_{w_{ij}}$ is solved as

[TABLE]

where $a={\frac{k_{e}q_{i}q_{j}}{2Td}}$ and ${b}={\frac{k_{e}|q_{i}q_{j}|}{2Td}}$. Putting $I_{w_{c}},I_{w_{ij}}$ in $I_{\mathbf{v}}^{res}$, and finally $I_{\mathbf{v}}^{res}$ in the restitutive part of Eq. (39), and integrating over $q_{i},q_{j}$, we obtain

[TABLE]

Notice that if $\mathcal{B}=0$, we recover the classical Haff’s law for viscoelastic granular gas.

IX Solution for the aggregative part $\left(\frac{3}{2}\frac{\partial T}{\partial t}\right)_{agg}$

The solution for the parts $I_{\mathbf{n}}^{agg},I_{\mathbf{v}}^{agg}$ are as follows.

[TABLE]

Using Eq. (52) and (36) from the above text, and again using the variable transformations (i) $\mathbf{w}_{ij}=\frac{2T}{m}(\mathbf{v}_{i}-\mathbf{v}_{j})$, (ii) $\mathbf{w}_{c}=\frac{2T}{m}(\mathbf{v}_{i}+\mathbf{v}_{j})$, (iii) $d\mathbf{v}_{i}d\mathbf{v}_{j}=-\frac{1}{8}\left(\frac{m}{4T}\right)^{-3}d\mathbf{w}_{ij}d\mathbf{w}_{c}$, the integral $I_{\mathbf{v}}^{agg}$ in (41) reduces to

[TABLE]

Here notice that now the relative velocity limits are from $v_{ij}=0$ to $\sqrt{b}$, the condition for aggregation selected by $\Theta(\sqrt{b}-v_{ij})$. Finally putting (53) into aggregative part of (39) and integrating over $q_{i},q_{j}$, we obtain

[TABLE]

Notice that after integrating from $v_{ij}=0$ to $\sqrt{b}$, the integration over $q_{i},q_{j}$ is to be broken into the sum of two parts, one over $q_{i}\in(-\infty,0]$, $q_{j}\in[0,+\infty)$, plus a second integral over $q_{i}\in[0,+\infty)$, $q_{j}\in(-\infty,0]$, to satisfy the aggregative condition set by $\Theta(-q_{i}q_{j})\Theta(\sqrt{b}-v_{ij})$.

A similar procedure is followed for $\left(\frac{\partial\langle\delta q^{2}\rangle}{\partial t}\right)_{res}$, $\left(\frac{\partial\langle\delta q^{2}\rangle}{\partial t}\right)_{agg}$ and $\left(\frac{\partial n}{\partial t}\right)_{agg}$ using corresponding $\Delta[\psi_{i}]$. Finally, the key constraint to be noted is that the solutions of the integrals in the rate equation for $T$ are real valued for $\mathcal{B}\leq 4$, while in equations for $\langle\delta q^{2}\rangle$ and $n$, they are real valued for $\mathcal{B}\leq 1$. The MD simulations confirm that these constraints put a physical limit during the aggregation phase.

X Granular MD simulations

The equation of motion of the form

[TABLE]

is solved for each particle with a setup of periodic boundary conditions in a cubic box of size $L^{3}=70d_{0}\times 70d_{0}\times 70d_{0}$, where $d_{0}$ is the monomer diameter [Table 2, 3]. Here $\mathbf{b}$ is a vector of integers representing the periodic replicas of the system in each Cartesian direction. The symbol ′ indicates that $k\neq i$ if $\mathbf{b}=0$ to avoid Coulomb interaction of particles with themselves. The non-dimensional numbers in the above equation are $\mathcal{E}=\frac{\alpha l_{\text{ref}}^{3/2}t_{\text{ref}}}{m_{\text{ref}}v_{\text{ref}}}$, $\mathcal{D}=\frac{\beta l_{\text{ref}}^{1/2}t_{\text{ref}}}{m_{\text{ref}}}$ and $\mathcal{K}=\frac{k_{e}q_{\text{ref}}^{2}t_{\text{ref}}}{m_{\text{ref}}v_{\text{ref}}l_{\text{ref}}^{2}}$, with $\alpha$ and $\beta$ being viscoelatic material constants. From practical problems, we select the reference length $l_{\textrm{ref}}\equiv d_{0}$, time reference $t_{\textrm{ref}}$, velocity reference $v_{\textrm{ref}}$, and charge reference $q_{\textrm{ref}}$ such that the elastic force strength $\mathcal{E}\approx 278$, dissipative force strength $\mathcal{E}/10$, and Coulomb force strength is varied across $\mathcal{K}=0.4-5.0$ [Table 2, 3]. The effect of dissipation compared to elastic forces is extensively studied for neutral systems Brilliantov and Pöschel (2010). The variation of Coulomb strength compared to dissipation and elastic forces we have repoted in Singh and Mazza (2018). The above order of magnitudes of $\mathcal{E}$, $\mathcal{D}$, and $\mathcal{K}$ also helps to attain an early clustering in non-dimensional time units in a finite size ($N\sim 50000$) neutral granular gas system [see Singh and Mazza (2018) for more details]. Also $\xi_{ij}\equiv d_{0}-r_{ij},\;\dot{\xi}_{ij}\equiv\frac{d\xi}{dt}$, $\mathbf{n}_{ji}$ is the unit vector pointing from center of in-contact neighbor $j$ towards the center of particle $i$, while $\mathbf{r}_{ki}$ is the distance vector pointing from particle $k$ towards the center of particle $i$.

The equation of charge on particle $i$ may be written as

[TABLE]

with $I_{ji}$ 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 step $\tau$ by Eq. (27), i.e.

[TABLE]

After charge-exchange, the long-range Coulomb forces for a setup with periodic boundary conditions in Eq. (X) 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 $\mathcal{O}(N^{3/2})$. 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 $\sim 10^{3}$ for a typical simulation with monomers $N\sim 10^{5}$, including the long-range electrostatic forces, is of the order of weeks. See Singh and Mazza (2018) for more details.

XI Comparison of individual $T$, $\langle\delta q^{2}\rangle$, and $n$ profiles, and emergence of convective flow using Mach number

In Fig. 7, we decompose the theoretical comparison of $\mathcal{B}$, presented in the main text, into individual comparisons of $T$, $\langle\delta q^{2}\rangle$, and $n$ profiles for a typical simulation run. The difference between the kinetic theory with and without aggregation is also emphasized.

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 $T$ 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 $\langle\delta q^{2}\rangle$ 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 event caused by the emergence of macroscopic flow. The theory predicts the decay of cluster density only in an average sense [Fig. 7(c) inset]. To quantify the emergence of the macroscopic flow, we calculate the local Mach number

[TABLE]

where $\mathbf{V}$ is the macroscopic velocity, and $v_{T}$ the thermal velocity. 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 ($tv_{\textrm{ref}}/d_{0}\approx 10^{2}$), indicating an intricate aggregation and fragmentation dynamics, and the generation of macroscopic flow.

XII REFERENCE SCALES AND Laboratory relevance of present results

Typically non-Brownian growth in planetary dust becomes dominant for monomer sizes near or above several $\mu\mathrm{m}$ Zsom et al. (2010) and the growth barrier problem Spahn and Seiß (2015) starts to arise near $d\sim 10^{-3}\;$\mathrm{m}. The mass of silica particles in this range of sizes is $m\sim 10^{-4}-10^{-6}\;$\mathrm{k}\mathrm{g}. If the particles are initially agitated with velocities $v\sim 1.0\;$\mathrm{m};\mathrm{s}^{-1}, the time scale reference to convert our simulation time to laboratory time is $d/v_{\mathrm{ref}}\sim 4.78\times 10^{-3}\;$\mathrm{s}. Thus in our results the growth over $10^{4}$ units of non-dimensional time approximately implies growth over $\sim 10\;$\mathrm{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 $\approx 2\;$\mathrm{m}\mathrm{m} to $\approx 7\;$\mathrm{c}\mathrm{m} in $\approx 10\;$\mathrm{s} time for particles of such size and mass, and for initial monomer filling fraction of $\phi=0.076$) [Fig. 8].

XIII Codes and media

The following computer codes and media are made public:

MATLAB program to solve Eq. (5)-(7). 2. 2.

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

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/cluster-analysis. 4. 4.

Movie of the evolving and aggregating structures in the charged granular gas.

XIV Data availability

All data are available from the authors.

XV 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.