© The Authors 2025

1 Introduction

Undeniably, a process exists that leads to the formation of planets. The common conception is that this process starts with the agglomeration of micrometre-sized dust grains within protoplanetary disks, which stick because of cohesive contact forces. However, many studies came to the conclusion that growth should not continue beyond sizes of millimetres to metres, based on our current understanding of agglomerate stability and dust migration. For example, the bouncing barrier is a consequence of increasingly elastic collisions at lower collision velocities (see Dominik & Tielens 1997; Tanaka et al. 2012; Zsom et al. 2010 and a comprehensive review by Birnstiel et al. 2016), while the metre-size barrier limits growth because of rapid inward drift and the depletion of solids (Weidenschilling 1977; Nakagawa et al. 1986). In addition to these theoretical and numerical studies, the bouncing barrier was also observed in a number of microgravity experiments (Blum & Wurm 2008; Weidling et al. 2009; Güttler et al. 2010; Kelling et al. 2014; Demirci et al. 2017).

Mechanisms that help overcome growth barriers have been discussed, such as streaming instabilities and pressure bumps (e.g. Youdin & Goodman 2005; Simon et al. 2016; Squire & Hopkins 2018; Pinilla et al. 2012), but additional processes may be at play. One such mechanism that is suggested to aid growth past the bouncing barrier is the attraction of grains due to electrical charging, mainly through the triboelectric effect (contact electrification). Recent studies have found that charged particles compose significantly more stable agglomerates (Ruan & Li 2022) and increase the likelihood of larger agglomerates (Steinpilz et al. 2020; Singh & Mazza 2018). Notably, even though Steinpilz et al. (2020) and Singh & Mazza (2018) studied the agglomerate size distribution (ASD) in different setups, both recovered a power-law distribution with approximately the same exponent.

To find the crucial mechanisms in early-stage planet formation within systems similar to protoplanetary disks, that is, very dilute and driven granular gases, a new approach to simulations is necessary. Previous methods can be roughly separated into two categories, which suffer from different shortcomings:

First, there are probabilistic simulations as performed by, for example, Birnstiel et al. (2010), Johansen et al. (2008), Estrada et al. (2016), and Zsom et al. (2010). These simulations typically feature a very detailed description of the disk parameters that influence agglomerate growth. However, the collisions themselves, an integral part of the growth process, are only described by the purely probabilistic Smoluchowski equation. Specifically, coagulation and fragmentation behaviour is governed by the kernel functions used in these simulations. Thus, the choice of kernels has a strong impact on the results, and any kernel of reasonable complexity can only represent the actual collision in an extremely simplified way – even before one tries to include electrostatic interactions, as done by, for example, Ivlev et al. (2002) and Okuzumi et al. (2011).

Second, in complete contrast, there are discrete element method (DEM) simulations of granular gases, for example, Steinpilz et al. (2020) and Singh & Mazza (2018). While agglomerate-agglomerate collisions are calculated in detail, it is practically impossible, or at least extremely inefficient, to implement an environment that resembles a protoplanetary disk. To understand why this is the case, consider the following: at typical particle densities, the mean free path of agglomerates is very long; consequently, the ASD evolution takes place on timescales longer than years. On the other hand, particle–particle interactions take place on timescales of nanoseconds. Therefore, the overwhelming majority of the simulation time is occupied by dynamics that do not need to be resolved by a DEM code; only for the collisions themselves is it mandatory. Putting the above into numbers, at a reasonable time step size of 10⁻⁶ s, simulating the evolution of just one year would take roughly 10¹³ time steps – an impossible amount considering that each step will include expensive computations. Another important shortcoming of the previous DEM studies is that only freely cooling systems were considered. The granular gas in a protoplanetary disk, however, is driven by a number of mechanisms, for example hydrodynamics and gravity.

Based on these considerations, we propose a novel hybrid method for ASD evolution simulations in which we use a kinetic Monte Carlo approach for the time evolution between two collisions. This is the natural probabilistic partner to the deterministic DEM because together they describe a specific trajectory across the disparate timescales. Even though kinetic Monte Carlo simulations could describe multi-agglomerate collisions, for a dilute granular gas it is justified to only consider binary collisions. Furthermore, since all other particles are very far away, the influence of any long-range interaction is negligible compared with the interactions between the collision partners themselves. Additionally, the coupling to an external driving mechanism, in the present case the coupling to the turbulent gas, is assumed to be strong enough to allow us to neglect long-range interactions for almost the whole time between collisions, removing the need to take into account corresponding correlations between agglomerates. Therefore, the hybrid method requires the granular gas to be sufficiently diluted such that the time between collisions exceeds the characteristic coupling time.

These assumptions allow for iterative simulations of binary collisions, which have low computational costs as the particle number and simulated times are relatively low and short, respectively. A code that is able to reproduce the results presented in Sect. 3 is publicly available¹.

2 Methods

The hybrid method described in this work focuses on the agglomeration behaviour of charged particles; the general concept is applicable to all kinds of interactions. In fact, only considering short-range interactions (e.g. adhesion) would simplify the method. The following description should be viewed as a blueprint; the details (interactions, driving mechanisms, modifications to collision rates etc.) may be altered to suit the needs for a plethora of disk models. This alteration would be simple, compared with changing a regular DEM simulation, as many processes only influence the setup for a collision, but not the actual particle dynamics during a collision. For example, including sedimentation as a driving mechanism would only change the formula used to calculate the relative velocities, discussed in Sect. 2.1.4.

2.1 Reducing granular gas dynamics to binary collisions

Using this method, we simulated the evolution of a population of agglomerates in a representative volume element of some region in the protoplanetary disk. From this population, we randomly chose two agglomerates for collision. The agglomerates resulting from the detailed simulation of this collision, whose numbers can possibly range anywhere from one (perfect merger of the collision partners) to the number of particles making up the collision partners (complete pulverisation), replace their parent agglomerates in the population. Then, the next agglomerate pair is chosen.

Each agglomerate consists of particles at certain coordinates with velocities, angular velocities (all in the agglomerate’s centre of mass coordinate system) and fixed charges in their centres. Hence, even two agglomerates of the same mass may drastically differ in porosity, shape, total charge etc.

As agglomerates can stick together in a collision, the number of agglomerates within the population may decrease. To prevent finite size effects, the population is self-duplicated (and the representative volume element accordingly becomes twice as big) every time it reaches half of its original size; therefore, also a population of half the original size must be a representative sample of the agglomerate distribution in the system we are looking to describe. This method was introduced by Liffman (1992, termed ‘topping up’) and further developed into the constant-N algorithm by Smith & Matsoukas (1998), which has been used to study agglomerate size evolution in protoplanetary disks (e.g. Ormel et al. 2007). In this work, however, we used Liffman’s topping up algorithm, as we do not know how sensitive the size evolution is with regard to changes to the size distribution (the fact that topping up does not alter the shape of the distribution is a key difference compared with the constant-N algorithm).

It should be noted that we use the term ‘agglomerate’ to also include monomers, special agglomerates that only consist of a single particle. When we use the term ‘particle’, however, we are referring to a primary particle, irrespective of what kind of agglomerate it might be a part of.

2.1.1 Macroscopic time

By the nature of this hybrid simulation technique, the time evolution of all results is given through the number of collisions that have occurred until a certain point. In general kinetic Monte Carlo terms, the collisions are events with given rates of occurrence; the overall theory is applicable to all kinds of such events (see e.g. Young & Elcock 1966 and Jansen 2012).

To convert the number of collisions into a real, physical time, each realised collision has to be associated with a waiting time δt, that is, the time between the previous and the current collision. It is distributed according to an exponential probability density function, $$f(\delta t)=\frac{1}{⟨\delta t⟩}\exp \left(\frac{\delta t}{⟨\delta t⟩}\right),$$(1)

where the average ⟨δt⟩ is the inverse of the sum of the rates for all possible collisions in the current state of the population. Assuming homogeneity within the simulated representative volume element V, the number density of each individual agglomerate is 1/V. Therefore, given the relative velocity v_ij and the collision cross-section σ_col,ij, we write the rate for the collision of an agglomerate i with a specific other agglomerate j as $${\mathrm{\Omega }}_{ij}=\frac{{\sigma }_{\mathrm{c}\mathrm{o}\mathrm{l},ij}{v}_{ij}}{V}.$$(2)

How big is the representative volume element for a given mass of solids m_s? In our simulation, all N particles have the same mass, m, so that m_s = N m. Within the scope of this paper, we assumed a fixed mass ratio of gas to solids β := m_g/m_s = 100 (e.g. Draine et al. 2007). Let ϱ_g = m_g/V_g be the gas density, where we can neglect the volume of the solids to write V_g ≈ V. Thus, we find the volume, V, in dependence of the total mass of solids, $$V=\frac{\beta m_{\mathrm{s}}}{{\varrho }_{\mathrm{g}}},$$(3)

and, consequently, the average waiting time as $$⟨\delta t⟩=\frac{V}{{\displaystyle \sum _{i>j}}{\sigma }_{\mathrm{c}\mathrm{o}\mathrm{l},ij}{v}_{ij}}=\beta \frac{m}{{\varrho }_{\mathrm{g}}}\frac{N}{{\displaystyle \sum _{i>j}}{\sigma }_{\mathrm{c}\mathrm{o}\mathrm{l},ij}{v}_{ij}},$$(4)

where the summation takes all possible collisions into account exactly once. As these rates generally change (and are recomputed) after each collision, the mean of the waiting time distribution f(δt) is not constant. Following this discussion, we can now identify the time of a specific collision as the sum of all previously sampled waiting times, δt.

2.1.2 Collision probability

Based on the results from the previous section, $$P_{ij}={\mathrm{\Omega }}_{ij}⟨\delta t⟩$$(5)

is the probability that agglomerates i and j collide in the time span ⟨δt⟩, given that a single collision occurs (we assumed that the collision probabilities are so small that higher orders [multiple collisions] are unimportant).

As Coulomb interactions between different agglomerates are negligible compared with the coupling to the gas for the majority of the time between two collisions, we assume that they have no significant impact on the relative velocity in the collision rate – an assumption that does not hold for the collision cross-section. Neglecting higher order terms, the collision cross-section of two charged spheres with radii ρ₁ and ρ₂, a reduced mass, $$M_{\mathrm{r}\mathrm{e}\mathrm{d}}=\frac{M_1{M}_2}{M_1+M_2},$$(6)

and charges Q₁ and Q₂ is given independently by Ivlev et al. (2002) and Scheffler & Wolf (2002) as $${\sigma }_{\mathrm{c}\mathrm{o}\mathrm{l}}=\pi {({\rho }_1+{\rho }_2)}^2(1-\frac{Q_1{Q}_2}{{\rho }_1+{\rho }_2}\frac{2}{M_{\mathrm{r}\mathrm{e}\mathrm{d}}{v}^2}).$$(7)

For non-spherical agglomerates, the radii ρ₁ and ρ₂ in the collision cross-section must be given by the maximal extent of the agglomerates, measured from the centre of mass. This is because, for anything less, the collision probability would not account for collisions with the maximum possible distance between the centres of mass. On the other hand, it is no problematic overestimation of the collision probability either, since it is still possible that the agglomerates do not collide in the DEM simulation. Therefore, this hybrid ansatz allows for the splitting of the collision probability into a trivial part (Eq. (7)) and a non-trivial part, which is directly resolved within the DEM simulation. The only drawback is that some unnecessary simulations are performed, that is, simulations in which the agglomerates do not come into contact, increasing the computation time. This is, however, a small price to pay with respect to the substantial benefits of this approach (i.e. not having to calculate all particle trajectories all the time for more than 10¹³ time steps).

2.1.3 Setup of imminent collisions

A number of quantities influence the initial setup for each collision. These are the initial distance between the agglomerates, their initial relative centre of mass velocity, the impact parameter, as well as their orientations.

The initial distance d_init determines the initial potential energy of the collision partners and, thus, has a direct relation to the impact velocities, as explained below. In the simulation, the distance is calculated as the maximal extent from the centre of mass of the larger agglomerate times a parameter, the initial distance coefficient, chosen to be 10 in the current work.

As long as d_init is chosen to be large enough that electrical dipoles and higher orders are negligible, the consequences of this choice can be absorbed into the initial relative velocity: for a (Coulomb interaction free) relative velocity v, the reduced mass m_red as well as the agglomerate charges Q₁ and Q₂, the initial velocity at distance d_init is given by $$v_{\text{init\hspace{0.17em}}}=\sqrt{v^2-2\frac{Q_1{Q}_2}{M_{\text{red}}{d}_{\text{init}}}}.$$(8)

In other words, the relative velocity can be split into a free relative velocity v, that is, the relative velocity in the absence of electrical charges, and the velocity change by converting potential energy from the distance travelled in the electric potential to kinetic energy.

Obviously, if the free relative velocity is too small to overcome the Coulomb barrier between two charged agglomerates of equal polarity, a collision is not possible at all. In this case, the respective collision probability is zero.

The simulation draws a random impact parameter b ∈ [0, b_max] for each collision from the probability density function $2b/b_{\text{max\hspace{0.17em}}}^2$, accounting for the fact that larger impact parameters are more likely in three dimensions. Using the simplest estimate, the simulation determines the maximum possible impact parameter from the collision cross-section as $b_{max}=\sqrt{{\sigma }_{\mathrm{c}\mathrm{o}\mathrm{l}}/\pi }.$.

The two agglomerates are initialised with random orientations. Their angular momenta follow from the particles’ positions and velocities. The initial population consists of dimers with random orientations and no angular momenta.

Even though all dynamical variables of the particles, including velocities, are kept between collisions, there will never be any significant amount of internal kinetic energy left from the last collision when two agglomerates make contact. This is because the dissipation timescales of the contact interactions are orders of magnitude shorter than d_init/v_init.

2.1.4 Relative velocities in a turbulent gas

An agglomerate of mass M couples via friction to the gas, which is present in a protoplanetary disk, with a characteristic stopping time $${\tau }_{\mathrm{s}}=\frac{3M}{4v_{\text{th}}{\varrho }_{\mathrm{g}}{\sigma }_{\text{aer}}}.$$(9)

Here, σ_aer is the aerodynamic cross-section, v_th the thermal gas velocity and ϱ_g the gas mass density (Ormel & Cuzzi 2007). The stopping time is related to the dimensionless Stokes number St = τ_s/t_L, where t_L is the orbital period, which we chose to be one year.

The aerodynamic (geometric) cross-section needed for the stopping time is a complicated quantity. Even assuming molecular free flow, which will break down for agglomerate sizes above centimetres (Weidenschilling & Cuzzi 1993), and instantaneous orientation of the agglomerates with respect to the gas flow, a unique orientation and, thus, a unique cross-section generally does not exist. To keep the approach simple, we used the radius of gyration to approximate the aerodynamic cross-section ${\sigma }_{\text{aer\hspace{0.17em}}}=\pi {\rho }_{\mathrm{g}\mathrm{y}\mathrm{r}}^2$. With the mass density ϱ_agg(r) of an agglomerate that has its centre of mass in the coordinate origin, the radius of gyration is given as $${\rho }_{\mathrm{g}\mathrm{y}\mathrm{r}}=\sqrt{\frac{1}{M}{\displaystyle \underset{{\mathbb{R}}^3}{\int }}{\mathrm{d}}^{3}r{\mathit{r}}^2{\varrho }_{\mathrm{a}\mathrm{g}\mathrm{g}}(\mathit{r})}.$$(10)

In the following, we discuss the model for turbulence-induced relative velocities, which greatly exceed the thermal, sedimentation and radial drift velocities (Blum et al. 1996). Here, the Reynolds number Re = 2.32 · 10⁷ and the velocity v_g = 18.5 m/s, which relates to the total kinetic energy of all eddies (see Eq. (2) in Ormel & Cuzzi 2007 and Eq. (46) in Birnstiel et al. 2011) characterise the turbulence. These values correspond to a turbulence strength of α = 10⁻⁴ (Weidenschilling & Cuzzi 1993).

The average squared relative velocity of two agglomerates (Δv)² := ⟨(v₁ – v₂)²⟩ is given by Ormel & Cuzzi (2007) as $$(\text{Δ}v{)}^2={\displaystyle \sum _{i\in \{1,2\}}}{v}_{\mathrm{g}}^2(1-G_i^{(1)}-G_i^{(2)}{G}_i^{(3)}),$$(11)

with the abbreviations $$G_i^{(1)}:=\frac{{\mathrm{S}\mathrm{t}}_i^2(1-{\mathrm{R}\mathrm{e}}^{-1/2})}{({\mathrm{S}\mathrm{t}}_i+1)({\mathrm{S}\mathrm{t}}_i+{\mathrm{R}\mathrm{e}}^{-1/2})},$$(12) $$G_i^{(2)}:=\frac{2{\tau }_{\mathrm{s},i}}{{\tau }_{\mathrm{s},1}+{\tau }_{\mathrm{s},2}},$$(13)

and $$\begin{array}{rl}{G}_i^{(3)}:=& 1-{\left(\frac{k^{\ast }}{k_{min}}\right)}^{-2/3}-\frac{{\mathrm{S}\mathrm{t}}_i}{1+{\mathrm{S}\mathrm{t}}_i}\\ & +\frac{{\mathrm{S}\mathrm{t}}_i}{1+{\mathrm{S}\mathrm{t}}_i{\left(k^{\ast }/k_{min}\right)}^{2/3}}.\end{array}$$(14)

Here, we use the index i ∈ {1, 2} to distinguish the values of the Stokes numbers and the stopping times of the two agglomerates considered as collision partners. The lowest wavenumber k_min is given by the inverse size of the largest eddy.

Furthermore, we need to discuss the wavenumber k*. It is given by $k^{\ast }=min\{k_1^{\ast },k_2^{\ast }\}$, where the wavenumber $k_i^{\ast }$ is the separation between two classes of eddies with respect to agglomerate i (see Völk et al. 1980):

Every agglomerate couples differently to the eddies, so this separation wavenumber also varies for different agglomerates. Importantly, the above formula, Eq. (11), is only defined for not too massive agglomerates, that is, only for agglomerates with τ_s < t_L, as the product $G_i^{(2)}{G}_i^{(3)}$ (the velocity correlation term) would be zero otherwise (k* = k_min). This is because all eddies could then be considered as class II eddies, which practically removes velocity correlations. Conversely, k*/k_min has the upper bound k_max/k_min = Re^3/4, for which all eddies are of class I. The highest wavenumber k_max is determined by the Kolmogorov scale. Ormel & Cuzzi (2007) give $$\frac{k_i^{\ast }}{k_{min}}=\frac{1}{2}{\mathrm{S}\mathrm{t}}_i^{-3/2}+1$$(15)

as an approximation for the separation wavenumbers.

The resulting relative velocities for typical stopping time values are depicted in Fig. 1. Notably, for a large relative velocity, at least one of the agglomerates needs to have a large stopping time; if both agglomerates are fully advected, they would only see class I eddies and would follow the streamlines perfectly.

Fig. 1 Relative velocities, Δv (see Eq. (11)), with respect to typical stopping times, τs,1 and τs,2, of the two collision partners.

2.2 Simulating a binary collision

To assume as little as possible about the results of a collision, we simulated it in the DEM framework, which is explained in the following.

2.2.1 Natural units

In this work we used a set of natural units that is based on Gaussian units. Notably, this means that the Coulomb force has no prefactor 1/4πε₀ and, as a consequence, that the unit of time can be expressed in units of length, mass, and charge, $$[\text{Time}]=\frac{\sqrt{[\text{Mass}][\text{Length}{]}^3}}{[\text{Charge}]}.$$(16)

Because of Eq. (16), many combinations of values whose units include at least three of the four above-mentioned ones are sufficient to construct a set of natural units. In this work specifically, we chose the particle diameter d, the particle mass m and a typical particle charge q_ref, which will serve as a reference charge. Thus, the corresponding natural time τ = 1 can be obtained from Eq. (16). In Table 1, the simulation parameters are listed. In order to give an idea of the orders of magnitude, plausible values in SI units are also included. We chose the particle diameter d as 0.1 mm, which is a typical size for particles used in experiments concerned with the growth and stability of agglomerates (e.g. Schwaak et al. 2024). The mass m is calculated from a particle mass density of 2500 kg/m³ (Blum et al. 1996). The reference charge of q_ref 10⁶ elementary charges is of the order of magnitude of observed tribocharges (e.g. Waitukaitis et al. 2014).

2.2.2 Contact force models

A number of contact interactions are implemented in the simulation to properly resolve the agglomerate collisions. These inter-actions only produce non-zero forces and torques when particles are in contact, that is, when the particles, all having diameter d, at positions r_i and r_j have an overlap ξ_n = d – ||r_i – r_j|| greater than zero. All interactions are assumed to be viscoelastic: included are a normal force – the Hertz–Kuwabara–Kono (HKK) model (see Hertz 1882; Kuwabara & Kono 1987; Brilliantov et al. 1996) and tangential, rolling and torsion frictions modelled by continuous spring-dashpot sliders (C-SDSs; see Führer et al. 2024). For a detailed description of the interactions, refer to Appendix A.

2.2.3 Electrostatic interactions

All particles carry a charge of magnitude q and sign s_i in their centre and interact via the Coulomb force: $${\mathit{F}}_{\mathrm{C}}=s_i{s}_j\frac{q^2}{{\parallel {\mathit{r}}_{ij}\parallel }^3}{\mathit{r}}_{ij}.$$(17)

Unfortunately, calculating all interactions between two large agglomerates is very costly, as the number of Coulomb forces grows quadratically with the number of particles; because of that, the complexity of a DEM time step for a collision with N_DEM particles is $O\left(N_{\text{DEM}}^2\right)$. Therefore, in the following, we discuss an approximation that reduces computational cost very efficiently.

Not describing an agglomerate as a number of primary particles but as a single rigid body with a given multipole field decreases the number of interactions to that of a two-body system – independent of the actual number of involved particles. For the simulations discussed in Sect. 3, we implemented the interactions for monopoles Q, dipoles p and quadrupoles Q. The formulas for the force F and the torque J on an agglomerate with multipole moments $\{\tilde{Q},\tilde{p},\tilde{\text{Q}}\}$ at the coordinate origin in the field of an agglomerate with multipole moments {Q, p, Q} at position r₀ are the result of a tedious but straightforward combination of the approximative electric field components $$E_{\nu }(\mathit{r})=\frac{Q{r}_{\nu }}{r^3}+3\frac{p_{\mu }{r}_{\mu }}{r^5}{r}_{\nu }-\frac{p_{\nu }}{r^3}+\frac{5}{2}\frac{r_{\mu }{Q}_{\mu \kappa }{r}_{\kappa }}{r^7}{r}_{\nu }-\frac{Q_{\nu \mu }{r}_{\mu }}{r^5}$$(18)

with the approximative force $$F_{\nu }=\tilde{Q}{E}_{\nu }(-{\mathit{r}}_0)+{{\tilde{p}}_{\mu }\frac{\mathrm{\partial }{E}_{\nu }}{\mathrm{\partial }{r}_{\mu }}|}_{-{\mathit{r}}_0}+{\frac{1}{6}{\tilde{Q}}_{\mu \kappa }\frac{{\mathrm{\partial }}^2{E}_{\nu }}{\mathrm{\partial }{r}_{\mu }\mathrm{\partial }{r}_{\kappa }}|}_{-{\mathit{r}}_0}$$(19)

or torque $$J_{\nu }={\epsilon }_{\mu \kappa \nu }{\tilde{p}}_{\mu }{E}_{\kappa }(-{\mathit{r}}_0)+{\frac{1}{3}{\epsilon }_{\mu \kappa \nu }{\tilde{Q}}_{\mu \lambda }\frac{\mathrm{\partial }{E}_{\kappa }}{\mathrm{\partial }{r}_{\lambda }}|}_{-{\mathit{r}}_0},$$(20)

respectively. Here, Greek letters as indices denote the Cartesian coordinates {x, y, z}, and monopoles Q cannot be confused with quadrupole tensor components Q_νμ because of the absence of indices. The full formulas are presented in Appendix B.

Three conditions govern the activation of the rigid-body approximation and ensure that this approximation is reasonable and efficient:

First, an agglomerate considered to be a rigid body is not allowed to come into contact with other agglomerates, as contact forces are not calculated for rigid bodies. Additionally, a description by a finite number of multipole moments is only an approximation of a charge distribution’s far field. To satisfy these points, a distance condition is implemented. It continuously checks whether a rigid agglomerate gets too close to other agglomerates, which leads to the unfreezing of the rigid agglomerate to a soft agglomerate, that is, a non-rigid object. Additionally, this condition also only allows soft agglomerates to freeze if they are sufficiently far away from all other agglomerates. The freezing-distance factor times the maximal extent of the larger of the agglomerates is the least distance between agglomerates that still allows for the rigid-body description: if two agglomerates get closer than this distance (and they are rigid), they unfreeze, and likewise, if they move further apart than this distance (and the other two conditions below are also fulfilled), they freeze. According to this definition, the critical distance is the same for freezing and unfreezing.

Second, even if an agglomerate satisfies the distance condition, it may be subject to internal movement such as vibrations or even restructuring. Simply freezing such an object might drastically change the time evolution of the system – which is why the temperature condition needs to be satisfied before a soft agglomerate may freeze. The name stems from the connection to the granular temperature of the agglomerate; however, the condition is stricter than simply demanding a small granular temperature. To calculate the internal movements, we changed into the body’s frame of reference, that is, we stripped the centre of mass motion and rotating with the agglomerate’s angular velocity. Then, the translational and rotational kinetic energies of each particle are calculated. The maximum particle energy – this is where the condition is more restrictive compared with a temperature threshold – is now monitored for a time period given by expected oscillation frequencies. Only if the maximum energy is lower than the freezing-temperature parameter for the entire time, an agglomerate is allowed to freeze.

Third, as the equations of motion as well as the multipole interactions of a rigid body are more complex than those of primary particles, it is to be expected that the benefits for freezing only outweigh the additional cost from a certain agglomerate size onwards. Therefore, the size condition only allows freezing for agglomerates that consist of at least as many particles as a chosen minimum freezing size. This sweet spot can be estimated by benchmarking the cost per time step. In contrast to the previous conditions, the size condition is purely performance related and has no further impact on the validity of the approximation.

The strictness of these conditions is controlled by three parameters: the freezing distance factor, the freezing temperature and the minimum freezing size. For the simulation results presented in this work, we chose a freezing distance factor and minimum freezing size of 3 and a freezing temperature of 10⁻³md²/τ². We checked on one of the datasets presented below that the parameters above yield the same results as a simulation with the rigid-body approximation turned off.

In addition to the performance increase through the reduction in the number of interactions, describing agglomerates as rigid bodies has another benefit. As the time step size is limited by the necessary resolution of particle collisions, as discussed in the next section, eliminating some or even all soft agglomerates may allow for larger time step sizes.

2.2.4 Time integration and step size

We used the midpoint method to integrate the DEM time evolution: First, we approximated all dynamical variables u_i(t) at time t + Δt/2 as $$u_i(t+\text{Δ}t/2)=u_i(t)+\frac{\text{Δ}t}{2}{\dot{u}}_i(t).$$(21)

Knowing these u_i(t + Δt/2), we either directly know their time derivatives ${\dot{u}}_i(t+\text{Δ}t/2)$ as well (e.g. velocities in the case of positions) or we can calculate them from the u_i(t + Δt/2) (e.g. forces and accelerations from positions and velocities). Then, the midpoint method yields $$u_i(t+\text{Δ}t)=u_i(t)+\text{Δ}t{\dot{u}}_i(t+\text{Δ}t/2)$$(22)

for the full time step t → t + Δt.

Importantly, each collision must be simulated long enough that no significant amount of re-agglomeration is omitted. As there is no reliable way to estimate whether agglomerates will collide again in a reasonable amount of time in many-particle systems (e.g. using escape velocities), we checked that our results do not change when simulating the collisions for twice as long, after all agglomerates start to pairwise move away from each other.

For the integration, the time step size Δt must not only be small to reduce the overall error, it is also important that it is smaller than the timescale of the fastest physical process that we want to simulate. In the case of this simulation, there are two important timescales: the contact time during a collision and the oscillation frequencies of the C-SDSs.

The contact time of two (uncharged) particles for given relative velocities is well described by a power law within the HKK model (see Schwager & Pöschel 1998). Therefore, its parameters depend on the stiffness k_n as well as the damping η_n and can be easily fitted from the results of a number of preliminarily simulated two-particle collisions at varying impact velocities. As the contact time decreases monotonically with increasing relative velocities, the shortest contact time is defined by the largest relative velocity between two particles that are in contact or may make contact in the next step. Given this minimum contact time τ_min, the current step size is set to the smaller one of $$\text{Δ}t=\frac{{\tau }_{\text{min}}}{20}$$(23)

and the maximum allowed step size Δt_max = 10⁻³τ, which must be small enough to properly resolve the oscillations of the C-SDSs. They behave like damped springs when oscillating; therefore, their frequencies only depend on the respective model parameters. The simulation uses Δt_max if all agglomerates are described as rigid bodies.

Fig. 2 Number of time steps for a collision simulation as a function of the total number of particles, NDEM, involved in the collision. Shown are 900 data points for collisions in the ensemble of populations with particle charge q = 1.1qref. Both axes have a logarithmic scale, and two power laws are shown for reference.

2.3 Computational cost and complexity

The simulations presented in Sect. 3 were performed on two Intel® Xeon® Gold 6348 CPU (2.60 GHz) processors, with run times of less than two weeks per simulated population. Each population’s evolution was calculated on a single core so that a total of 56 evolutions could be simulated at once.

Roughly 95% of the CPU time is spent on performing the DEM collisions, everything else (updating collision rates, drawing new collision partners, writing outputs etc.) only makes up the remaining few percent. Not using any further approximations, the complexity of a DEM time step involving N_DEM particles will be $O\left(N_{\text{DEM}}^2\right)$ for at least some time steps because of the long-range Coulomb interaction (the rigid body approximation will never be applicable during the actual contact of the agglomerates).

The number of DEM time steps necessary to carry out an agglomerate collision depends on, for example, the initial distance, the relative velocity, the agglomerates’ contact time, and the finishing conditions, all of which depend on the particle number in a complicated way. Hence, we present the number of time steps for some exemplary collisions in the investigated system in Fig. 2.

The complexity of a brute force update of the collision rates discussed in Sect. 2.1.1 for N_A agglomerates is $O\left(N_{\mathrm{A}}^2\right)$, but can be reduced to O (N_A) if the rates are only recalculated for collisions with agglomerates created in the most recent collision (assuming their number is of the order of one) for the system studied in this work. In general, if the environment (density, turbulence strength etc.) is not in a steady state, all collision rates have to be recalculated after each collision. The calculations of properties of an agglomerate of size M required for the rate update (maximal extent, radius of gyration etc.) have a complexity linear in M – however, as these quantities also have to be calculated in every time step of the DEM part, they are readily available and do not cost any additional time.

How to estimate the number of collisions needed to simulate a given physical time? The average waiting time, which is the average time increment per collision, is given by Eq. (4). The sum in the denominator has $N_{\mathrm{A}}^2,N_{\mathrm{A}}\in (N_{\mathrm{A},\text{init}}/2,N]$, terms (kept larger than N_A,init/2 by the topping up algorithm), where N_A,init denotes the initial number of agglomerates. This leads to a number of needed collisions of the order of $N_{\mathrm{A},\text{init}}^2$ or N² at best or worst, respectively. Moreover, the number of particles in the population N, appears in the numerator of Eq. (4), and also non-trivially affects the collision cross-sections and relative velocities in the denominator, as discussed earlier in this section. The situation is further complicated, as N and N_A change in time because of the topping up algorithm, agglomeration, and fragmentation. To get a quantitative idea about how many collisions were necessary for the results in Sect. 3, see Fig. 3. As expected, we notice a correlation between the number of collisions and the number of agglomerates in the populations at the final time t = 98 years. However, the influence of the number of particles is less clear. To understand Fig. 3 in more detail, we need some results from – and, thus, delay its discussion to the end of – Sect. 3.1.

Fig. 3 Number of collisions simulated and number of agglomerates at the end of the evolution at t = 98 years for each of the 3 × 56 populations discussed in Sect. 3. All populations started with NA,init = 500 dimers. The colours indicate the three ensembles, which differ in terms of particle charges, q. The populations with q = 0.9qref finished with N = 2000 particles (duplicated once) and the populations with q = 1.0qref finished with N = 4000 particles (duplicated twice). Two distinct sets of populations with q = 1.1qref are visible, with (NA < 325, N = 4000) and (NA > 475, N = 8000), respectively. The populations with N = 8000 are duplicated three times.

3 Results

In the presented simulations, the initial agglomerate population consists of neutral dimers with random orientations and no angular momenta. Each dimer is built out of two spherical particles in contact, which carry a charge of equal magnitude. Apart from their polarity, all particles are equal. We simulated an ensemble of 56 systems – differing in their (random) initial configuration – for each studied particle charge. Table 1 contains all other required parameters.

Each population started off with 1000 particles. To improve the statistics in the following analysis, we view the collection of all individual populations – which all are representative samples that evolved independently under equal conditions – as a single large population. As the macroscopic time step sizes (see Sect. 2.1.1) are stochastic in nature and interpolation is not an option for the given data, it may not be trivial to see, which individual populations are included in the combined population at a time t. Considering that we know the population exactly at the times, at which it changes (because of a collision), we also know that the population remains unchanged in the time between two collisions. Therefore, the combined population at time t consists of all those individual populations that were produced by the last collisions before the time t.

We want to point out that the term size distribution is often used for the distribution of the agglomerate masses (measured in numbers of constituent particles) and not for the distribution of, say, their radii (see e.g. Steinpilz et al. 2020). To be consistent with the afore mentioned literature, we too call the distribution of masses ASD.

Fig. 4 Time evolution of the average agglomerate size for three different particle charges, collapsed onto each other to reveal the universal exponential convergence. The average sizes, ⟨M⟩(t), and average asymptotic sizes, ⟨M⟩∞, are measured in numbers of primary particles. Both the average asymptotic sizes and the characteristic times, t*, are used as fit parameters; refer to Table 2 for their values. The dashed line is included to guide the eye.

3.1 Stationarity of the size distribution

As we simulate a driven granular gas where agglomeration and fragmentation of agglomerates can occur, the possibility of an emerging stationary ASD is a reasonable assumption. For example, Zsom et al. (2010) discussed the possibility of emerging quasi-steady states in the midplane of a disk due to bouncing, and Brilliantov et al. (2015) explained the stationarity of the grain size distribution in Saturn’s rings by assuming an agglomeration–fragmentation equilibrium. Here we should distinguish between the two trivial (and extreme) distributions – the existence of either only monomers or only one agglomerate that contains the entire mass of the population – and all other non-trivial ones. In the context of planet formation, conditions that result in the latter trivial case, that is, conditions that result in a run-away growth, are of particular interest.

To check whether the ASD reached a stationary state, it is, in principle, necessary to check if all moments of the distribution become stationary. In practice, however, this is not feasible. For the purpose of this work, we checked the evolution of the average agglomerate size (Fig. 4) and a contour plot of the entire distribution’s evolution (Fig. 5). Both figures indicate stationarity of the ASD for the two larger charges; however, the simulated time was not sufficiently long to reach the stationary state of the system with q = 0.9q_ref, which converges more slowly.

Greater charges result in larger agglomerates – as might be expected – visible in Fig. 5 and quantified by the asymptotic average agglomerate sizes ⟨M⟩_∞ in Table 2. In addition, the time needed to reach the stationary distribution depends on the particle charge as well: the characteristic time t* of the average size’s exponential convergence to its asymptotic value gets smaller with increased charge (see Table 2).

Now that we know that the average agglomerate size, which is given by the number of particles, N (equivalent to the total mass of the population in our case) and the number of agglomerates, N_A, by definition, $$⟨M⟩=\frac{N}{N_{\mathrm{A}}}$$(24)

converges to an asymptotic value ⟨M⟩_∞ (see Fig. 4), we can understand Fig. 3 in more detail. Using our fit values for ⟨M⟩_∞ (Table 2), we can calculate asymptotic values for the ensemble average of the number of agglomerates N_A,∞(N). A single population of finite size will always only fluctuate around this value, and only the ensemble average will actually tend to it. For example, we get N_A,∞ ≈ 272 and N_A,∞ ≈ 544 for N = 4000 and N = 8000, respectively, for the populations with particle charges q = 1.1q_ref. This is consistent with with the data shown in Fig. 3: the populations with particle charges q = 1.1q_ref are split into two sets, according to their numbers of particles N = 4000, 8000. In general, our populations start off with very small agglomerates, which leads to a decreasing number of agglomerates due to agglomeration. Neglecting fluctuations, the topping up algorithm will keep duplicating the populations, until N_A,∞(N) is larger than the lower bound for the number of agglomerates, N_A,init/2 = 250. In the case of the q = 1.1q_ref-ensemble, however, the asymptotic number of agglomerates N_A,∞(4000) ≈ 272 is so close to the lower bound that fluctuations have caused some populations to be duplicated, resulting in the two sets of populations visible in Fig. 3.

Fig. 5 Time evolution of the ASD for different particle charges. The agglomerate sizes are measured in numbers of primary particles, and their frequencies are normalised by the total number of agglomerates at a certain time.

3.2 Agglomerate size distribution

As previously stated, both Steinpilz et al. (2020) and Singh & Mazza (2018) studied the ASD for charged particles – in contrast to this work, for freely cooling systems – and both reported a power law with an exponent of −1.3 to −1.4, which roughly fits with the exponent of −1.5 produced in our simulations (see Fig. 6). However, the power law describes the ASD only for small agglomerates: the tails of the distribution actually follow an exponential distribution. This is also supported by the fact that the deviation from the power law can already be seen in the results from Singh & Mazza (2018) (in their Fig. 4, the data deviate from the power-law fit for s > 10; s is their measure for the agglomerate size). We estimated the probabilities by calculating the difference quotients of the empirical cumulative distribution function for agglomerate sizes. The employed difference quotients only take into account sizes that occur at least once.

By rescaling the agglomerate sizes, we collapsed the ASDs in Fig. 6 onto each other and show that the distribution changes under the variation of the particle charge only through the characteristic size M*, which we used as a fit parameter (see Table 2). Notably, this characteristic size M* is larger than the average size ⟨M⟩(t = 98 yrs) in Fig. 4. This is because small agglomerates (the part that fits the power-law description) are more frequent than a pure exponential distribution would suggest. Therefore, the average is skewed towards smaller agglomerates, whereas the exponential decay in the tails is not affected.

4 Conclusion

In this paper we have presented a novel hybrid method for simulations of ASD evolutions, combining the efficiency of probability-based simulations and the precision of DEM simulations for the microscopic processes of agglomerate collisions. We find that the average agglomerate size converges exponentially fast to a charge-dependent asymptotic size, where the increased stability of higher charges allows for further agglomeration. Furthermore, an analysis of the full ASD reveals that the system reaches a steady state on a timescale of 100 years for our parameter choices.

The observed exponential tails of the ASD are seemingly in contradiction with the power laws presented by Steinpilz et al. (2020) and Singh & Mazza (2018). However, the present simulation was able to produce better statistics for agglomerates of far larger particle numbers, which allowed us to study the tail behaviour in much more detail. Our findings are supported by the visible deviation from the Singh & Mazza (2018) power law shown in their Fig. 4, panels b–d. From these results alone, it is not possible to make any quantitative claims about size distributions in protoplanetary disks; more extensive parameter studies are needed. However, the presence of an exponential drop-off in the steady state might be an indication of the existence of a growth barrier, beyond which agglomerate sizes are not impossible but extremely unlikely.

We want to point out that it generally is straightforward to extend the presented model to include additional physical processes. As long as the main assumptions – a sufficiently dilute granular gas and collision times much shorter than particleenvironment coupling times – are satisfied, the effects of different driving mechanisms (e.g. radial drift or sedimentation) and spatial inhomogeneities of the agglomerate number density, just to name a few, solely affect the relative velocities and collision rates. Therefore, a slight modification of the code allows for the study of those processes’ roles in a plethora of systems.

While in this paper we have focused on the role of the particle charges, studying local variations in turbulence strength might be crucial in understanding the local disk structure with regard to the ASD. Furthermore, an investigation of the robustness of the qualitative behaviour and the shape of the ASD under different driving mechanisms would be an important continuation of this study. Also, it is possible to classify collisions according to their outcomes (e.g. the types proposed by Zsom et al. (2010), their Sect. 3.4.) and use these statistics to evaluate their relevance with respect to agglomerate growth and growth barriers. Additionally, as the results of our hybrid method encompass not just agglomerate sizes but also the actual positions of the constituent particles, we intend to analyse these collisionally grown agglomerates with respect to porosities and binding energies for charged as well as adhesive particles (long-range vs short-range interactions) in the future.

Fig. 6 ASDs for three different charges at t = 98 years. All agglomerate sizes, M, are rescaled with characteristic sizes, M* (see Table 2), such that the tails have a slope of −1 in the semi-logarithmic (lower) plot. The exponential function and power law are included to guide the eye. The upper and lower plots show the same data but with logarithmically and linearly scaled abscissae, respectively.

Acknowledgements

We thank the anonymous referee for thoroughly examining our manuscript and for suggesting many changes that strengthened the message of this paper. The presented research is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer(458889524). We acknowledge support by the Open Access Publication Fund of the University of Duisburg-Essen.

Appendix A Details of implemented contact forces

All contact forces have in common that they only act on particles that are in direct contact, that is, with a particle radius R and particle positions r_i and r_j, the overlap ξ_n = 2R – ||r_i – r_j|| must be greater than zero. For this study, all particles are of the same size, 2R, and mass, m.

We note that we did not implement any cohesive forces. This is because, as already mentioned by, for example, Steinpilz et al. (2020), the cohesion between micrometre-sized particles is way too small to have any significant influence on the stability of agglomerates when collision velocities exceed cm/s. Therefore, the cohesion would always be negligible compared with the Coulomb forces presented in Sect. 2.2.3.

We denote the spring stiffness and damping parameters as well as friction coefficients as k, γ and μ, respectively. We add an index to these parameters, indicating to which type of contact force they belong (see Table 1): n (normal), t (tangential), s (torsion or spin), and r (rolling). We chose all parameters such that the contact is critically damped, and the resulting time steps yield reasonable computational efforts.

A.1 Normal force

The normal force presented in the following is denoted as the force acting on the particle i; thus, the normal vector is $$\mathit{n}=\frac{{\mathit{r}}_i-{\mathit{r}}_j}{\parallel {\mathit{r}}_i-{\mathit{r}}_j\parallel }.$$(A.1)

Furthermore, we used the effective radius $$R_{\text{eff\hspace{0.17em}}}=\frac{RR}{R+R}=\frac{R}{2}$$(A.2)

and the analogously defined effective mass $$m_{\text{eff\hspace{0.17em}}}=\frac{m}{2}.$$(A.3)

The relative velocity v_ij = v_i – v_j has the normal component $${\mathit{v}}_{\mathrm{n}}=v_{\mathrm{n}}\mathit{n}=(\mathit{n}\cdot {\mathit{v}}_{ij})\mathit{n}.$$(A.4)

For the normal contact force F_n, we implemented the HKK model (Hertz 1882; Kuwabara & Kono 1987; Brilliantov et al. 1996): The repulsion term (Hertz) and the damping term (Kuwabara–Kono) are combined into the total normal contact force $$F_{\mathrm{n}}=\sqrt{R_{\mathrm{e}\mathrm{f}\mathrm{f}}}{k}_{\mathrm{n}}{\xi }_{\mathrm{n}}^{3/2}-m_{\mathrm{e}\mathrm{f}\mathrm{f}}\sqrt{R_{\mathrm{e}\mathrm{f}\mathrm{f}}{\xi }_{\mathrm{n}}}{\gamma }_{\mathrm{n}}{v}_{\mathrm{n}}.$$(A.5)

A.2 Tangential forces and torques

All implemented tangential forces and torques include sticking and sliding, implemented via the continuous spring-dashpot slider (Führer et al. 2024), which is a modification of the discrete spring-dashpot slider as presented by Luding (2008). Apart from the corrected spring-displacement update, the forces, torques, and relative velocities are implemented analogously to in Luding (2008). In essence, all forces and torques are calculated from the limited reaction force of a deforming viscoelastic material, $${\mathit{F}}_{\mathrm{S}\mathrm{D}\mathrm{S}}=min(-k\mathit{\xi }-\gamma \mathit{v},\mu F_{\mathrm{n}}),$$(A.6)

where ξ denotes the spring displacement, modelling the material deformation, and v the respective relative velocity, not considering material deformation (for details, see Führer et al. 2024).

We note that the parameters k, γ, and μ differ in their units for the different modes; only tangential and rolling interactions share the same units. This is because Hertz’s contact model is based on the contact area and not a linear displacement, and the torsion mode is concerned with angular (not length) changes.

Appendix B Multipole interaction formulas

The full formulas for the second order multipole–multipole interactions, calculated from Eqs. (18), (19) and (20) — using $$A=Q_{v\mu }{\tilde{Q}}_{v\mu },$$(B.1) $$(\mathit{B}{)}_{\kappa }={\epsilon }_{\nu \mu \kappa }{\tilde{Q}}_{\nu \mu },$$(B.2) $$(\mathit{C}{)}_{\kappa }={\epsilon }_{\nu \mu \kappa }{\tilde{Q}}_{\nu \lambda }{Q}_{\lambda \mu },$$(B.3)

and the identity matrix I — are $$\begin{array}{rl}& \mathit{F}=\frac{\tilde{Q}}{r^3}(Q\mathit{r}+3\frac{\mathit{p}\cdot \mathit{r}}{r^2}\mathit{r}-\mathit{p}+\frac{5}{2}\frac{\mathit{r}\cdot \mathrm{Q}\mathit{r}}{r^4}\mathit{r}-\frac{\mathrm{Q}\mathit{r}}{r^2})\\ & +\frac{Q}{r^3}(1-3\frac{\mathit{r}\otimes \mathit{r}}{r^2})\tilde{\mathit{p}}+3\frac{1}{r^5}(\mathit{p}\cdot \tilde{\mathit{p}}\mathit{r}+\mathit{p}\cdot \mathit{r}\tilde{\mathit{p}}-5\frac{\mathit{p}\cdot \mathit{r}\tilde{\mathit{p}}\cdot \mathit{r}\mathit{r}}{r^2})\\ & +3\frac{\tilde{p}\cdot \mathit{r}}{r^5}\mathit{p}+5\frac{\mathit{r}}{r^7}\cdot (\mathrm{Q}\tilde{\mathit{p}}\mathit{r}+\frac{1}{2}\mathrm{Q}r\tilde{p}-\frac{7}{2}\frac{\tilde{p}\cdot \mathit{r}}{r^2}\mathrm{Q}\mathit{r}\mathit{r})\\ & -(\mathrm{Q}\tilde{p}-5\frac{r\cdot \tilde{p}\mathrm{Q}\mathit{r}}{r^2})\frac{1}{r^5}+\frac{3}{6}\frac{Q}{r^5}(-2\tilde{\mathrm{Q}}\mathit{r}+5\mathit{r}\frac{\mathit{r}\cdot \tilde{\mathrm{Q}}\mathit{r}}{r^2})\\ & +\frac{3}{2r^5}(2\tilde{\mathrm{Q}}\mathit{p}-10\frac{\mathit{p}\cdot \tilde{\mathrm{Q}}r}{r^2}\mathit{R}+\mathit{p}\cdot \mathit{r}(-10\frac{\tilde{\mathrm{Q}}\mathit{r}}{r^2}+35\frac{\mathit{R}\cdot \tilde{\mathrm{Q}}r}{r^4}r)\\ & -5\frac{r\cdot \tilde{\mathrm{Q}}\mathit{r}}{r^2}\mathit{p})+\frac{5}{6r^7}(Ar+2\tilde{\mathrm{Q}}\mathrm{Q}r-14\frac{r\cdot \mathrm{Q}\mathrm{Q}r}{r^2}\mathit{r}\\ & +\frac{7}{2}r\cdot \mathrm{Q}r(-2\frac{\tilde{\mathrm{Q}}r}{r^2}+9\frac{r\cdot \tilde{\mathrm{Q}}r}{r^4}r)+2\mathrm{Q}\tilde{\mathrm{Q}}r-7\frac{r\cdot \tilde{\mathrm{Q}}r}{r^2}\mathrm{Q}r)\end{array}$$

and $$\begin{array}{rl}\mathit{J}=& \frac{\tilde{\mathit{p}}\times \mathit{r}}{r^3}(Q+3\frac{\mathit{p}\cdot \mathit{r}}{r^2}+\frac{5}{2}\frac{\mathit{r}\cdot \mathrm{Q}\mathit{r}}{r^4})-\frac{\tilde{\mathit{p}}\times \mathit{p}}{r^3}-\frac{\tilde{\mathit{p}}\times (\mathrm{Q}\mathit{r})}{r^5}\\ & +\frac{Q}{3r^3}(\mathit{B}-3\frac{(\tilde{\mathrm{Q}}\mathit{r})\times \mathit{r}}{r^2})+\frac{1}{r^5}((\tilde{\mathrm{Q}}\mathit{p})\times \mathit{r}+(\tilde{\mathrm{Q}}\mathit{r})\times \mathit{p}\\ & +\mathit{p}\cdot \mathit{r}(\mathit{B}-5\frac{(\tilde{\mathrm{Q}}\mathit{r})\times \mathit{r}}{r^2}))+\frac{5}{3r^7}((\widetilde{\mathrm{Q}\mathrm{Q}}\mathit{r})\times \mathit{r}+\frac{1}{2}\mathit{r}\cdot \mathrm{Q}\mathit{r}(\mathit{B}\\ & -7\frac{(\tilde{\mathrm{Q}}\mathit{r})\times \mathit{r}}{r^2})+(\tilde{\mathrm{Q}}\mathit{r})\times (\mathrm{Q}\mathit{r}))-\frac{\mathit{C}}{r^5},\end{array}$$

evaluated at r = −r₀.

References

All Tables

All Figures

	Fig. 1 Relative velocities, Δv (see Eq. (11)), with respect to typical stopping times, τs,1 and τs,2, of the two collision partners.
In the text

	Fig. 2 Number of time steps for a collision simulation as a function of the total number of particles, NDEM, involved in the collision. Shown are 900 data points for collisions in the ensemble of populations with particle charge q = 1.1qref. Both axes have a logarithmic scale, and two power laws are shown for reference.
In the text

Fig. 3 Number of collisions simulated and number of agglomerates at the end of the evolution at t = 98 years for each of the 3 × 56 populations discussed in Sect. 3. All populations started with NA,init = 500 dimers. The colours indicate the three ensembles, which differ in terms of particle charges, q. The populations with q = 0.9qref finished with N = 2000 particles (duplicated once) and the populations with q = 1.0qref finished with N = 4000 particles (duplicated twice). Two distinct sets of populations with q = 1.1qref are visible, with (NA < 325, N = 4000) and (NA > 475, N = 8000), respectively. The populations with N = 8000 are duplicated three times.

In the text

Fig. 4 Time evolution of the average agglomerate size for three different particle charges, collapsed onto each other to reveal the universal exponential convergence. The average sizes, ⟨M⟩(t), and average asymptotic sizes, ⟨M⟩∞, are measured in numbers of primary particles. Both the average asymptotic sizes and the characteristic times, t*, are used as fit parameters; refer to Table 2 for their values. The dashed line is included to guide the eye.

In the text

	Fig. 5 Time evolution of the ASD for different particle charges. The agglomerate sizes are measured in numbers of primary particles, and their frequencies are normalised by the total number of agglomerates at a certain time.
In the text

Fig. 6 ASDs for three different charges at t = 98 years. All agglomerate sizes, M, are rescaled with characteristic sizes, M* (see Table 2), such that the tails have a slope of −1 in the semi-logarithmic (lower) plot. The exponential function and power law are included to guide the eye. The upper and lower plots show the same data but with logarithmically and linearly scaled abscissae, respectively.

In the text