© The Authors 2023

1 Introduction

Understanding the collisional evolution of submillimeter-sized granular dust aggregates is crucial in many astrophysical environments. These aggregates are groupings of smaller individual grains (usually on micron scales) that are randomly interconnected. A grain (also called a monomer) is the smallest unit of study and can be assumed to be solid and homogeneous in composition. These grains are created via condensation in the protoplanetary disk of the Solar System; they were also previously created in the interstellar medium (ISM) or by asymptotic giant branch stars (Alexander et al. 2007).

Dust size plays a fundamental role in the interstellar processes (Mathis et al. 1977) that determine the state of the ISM, affecting its thermal and chemical equilibrium. According to Ormel et al. (2009), it is clear that the contents of molecular clouds, as well as the general state of the process of the formation of stars and planets, are linked to the properties of dust grains and, in particular, their size distribution. Recent space missions have provided relevant information on the properties of dust in dense clouds. In particular, comparisons of far-IR emission maps taken by the Infrared Astronomical Satellite (IRAS) and Spitzer satellite and the near-IR extinction maps derived from the Two Micron All-Sky Survey (2MASS) suggest grain growth in regions of higher density (Schnee et al. 2008). Coagulation processes have also been detected via a detailed absorption profile analysis of the 10 µm silicate absorption band in these environments (van Breemen et al. 2011).

In protoplanetary disks, collisions between dust aggregates can lead to the formation of larger aggregates, which are important intermediaries in planet building (Armitage 2010; Blum 2010). Even though dust growth does play a major role in shaping the evolution of protoplanetary dust and planet formation, it is sometimes neglected when building models of protoplan-etary disks due to its complexity and computational expense (Draczkowska et al. 2019). However, in recent years the upgrade of the existing (sub)millimeter arrays, the Atacama Large Millimeter/submillimeter Array (ALMA) and Very Large Array (VLA), has significantly improved the observational constraints on models of dust evolution in protoplanetary disks. Laboratory experiments and numerical simulations have led to a substantial improvement in the understanding of the physical processes of grain-grain collisions, which are the foundation for models of dust evolution in disks. Testi et al. (2014) review the constraints on the physics of grain-grain collisions as they have emerged from laboratory experiments and numerical computations; they also review the current theoretical understanding of the global processes governing the evolution of solids in protoplanetary disks, including dust settling, growth, and radial transport. At this moment, there is no single model that alone can explain all the aspects of the pathway from dust grains to planetary systems; Drazkowska et al. (2022) present an updated summary of planet formation theories, from dust evolution to the growth of planetary cores via the accretion of planetesimals, pebbles, and gas. Dust aggregates also make up the rings of planets and other minor objects (Burns et al. 2001). According to Burns et al. (1980), the small grains in Jupiter’s rings are generated by the impacts of micrometeorites, by collisions between particles of different sizes, and by self-fracture due to electrostatic stress.

Another important context is when these dust grains collide with larger objects, such as micrometeorites, that impact objects in the Solar System, including the Earth. This has allowed their study and classification (Genge et al. 2008). Also, these small particles can produce craters on the surfaces they impact and eject material from them. Likewise, when they hit a space probe or satellite, they can cause irreversible damage (Grossman et al. 2010). Finally, these types of collisions are also present in debris disks, which are remnants of the times of formation of planetary systems. Debris disks are believed to be created and maintained by mutual collisions, and (possibly) via cometary-like activity of some minor bodies, similar to asteroids, comets, and objects in the Solar System’s Kuiper belt. Dust evolves through a cascade of collisions under the action of stellar gravity and radiation forces (Krivov et al. 2006).

Making a complete dust model including growth and fracture (applicable to the scenarios described above) is not simple, and it should include the temporal evolution of the dust as its aggregates collide with one another, taking all possible parameters into account: different compositions, impact velocities, porosities, and size, among others. Attempts have been made to achieve this by including some of the most important parameters. However, the porosity of dust aggregates has rarely been taken into account in previous models because it is difficult to determine accurately since the scattering properties of large porous grains are qualitatively similar to those of small compact grains of the same mass and composition (Graham et al. 2007; Shen et al. 2009; Kirchschlager & Wolf 2013).

Gáspár et al. (2012) summarized the collision models applicable to debris disks and presented a code involving erosive and catastrophic collisions, employing an efficient numerical algorithm that allowed them to evaluate scattering integrals with high precision. Birnstiel et al. (2012b) derived a simple model for protoplanetary disks that describes the radial evolution of the dust surface density under the combined influence of the growth and fragmentation of compact grains as well as radial transport mechanisms. The important parameters are the fragmentation threshold velocity, the level of turbulence, the initial dust-to-gas ratio, and the temperature and density profile of the gas disk. They estimate the effective dust transport velocity by representing the dust distribution with only two characteristic grain sizes, a small and a large population. Drazkowska et al. (2013) developed a numerical model that resolves the spatial distribution of dust in the radial and vertical dimensions. The coagulation and fragmentation of solids are taken into account via a Monte Carlo method. A collision model that adopts the mass transfer effect (which can occur for different-sized dust aggregate collisions) is implemented. They focus on a proto-planetary disk that includes a pressure bump caused by a steep decline in turbulent viscosity around the snow line. None of the above three models involve aggregate porosity, but some authors have included the importance of porosity in the modeling of pro-toplanetary disks. Ormel et al. (2007) found that the evolution of their porous aggregates in a collisional model is quantitatively different from aggregation models in which porosity can be parameterized by a fixed exponent. In this context, coagulation models require a microphysical collision model that takes aggregate porosities into account. In a later work, Ormel et al. (2009) proposed a Monte Carlo code to study the collisional evolution of the dust population by combining two models: a binary model that simulates the collision between two aggregates and a coagulation model that computes the dust-size distribution with time. The collision model features the sticking, breakage, erosion, and shattering outcomes, and it includes off-center collisions and changes in the internal structure in terms of the porosity. It also allows for a scaling of the results to the relevant masses and critical energies, enabling the calculation to proceed beyond the sizes covered by the original numerical collision experiments. However, their binary model only includes collisions between aggregates with n ≤ 10³ particles, and they needed to extrapolate the outcomes for larger aggregates. Zsom & Dullemond (2008) also developed a Monte Carlo method that includes the porosity analysis of Ormel & Cuzzi (2007), in which they follow the history of a limited number of representative particles. Their code is almost ten times faster than the code by Ormel et al. (2007) because the calculation time scales linearly with the number of collisions simulated. The difference is fragmentation: in the simulations of Ormel et al. (2007), no fragmentation occurs because the growth timescales are longer, and the porosities of these particles would be smaller if fragmentation were included. This model can track only relatively few particles, and a limitation is encountered when strong growth and fragmentation occur simultaneously. Zsom et al. (2010) affirm that a more realistic collision model that includes the compaction and fragmentation of aggregates is necessary. They also used a Monte Carlo code to follow the mass and porosity evolution of the particles in time. They found that significant changes in the porosity of the aggregates have the potential to significantly alter their collision model and, therefore, the results obtained.

The evolution of the aggregate sizes and porosities has also been studied using coagulation models. Okuzumi et al. (2012) used the coupled equations of aggregate size and porosity and found that the growth of large, highly porous aggregates is possible if collisional fragmentation was negligible. Zsom et al. (2011) studied this combined evolution using a Monte Carlo scheme but pointed out that the evolution of porosity in particular required more detailed information (experimental or theoretical) to be reliable. Krijt et al. (2015) followed a similar strategy and concluded that the growth of porous particles showed considerable differences from the growth of compact particles. This work was later expanded to a semi-analytical model (Krijt et al. 2016). Homma & Nakamoto (2018) later pointed out that the radial drift in protoplanetary disks may play an important role in the assessment of growth and porosity evolution models, limiting the possibility of collisional growth. A coupling of radial drift to grain size and porosity evolution was implemented by Garcia & Gonzalez (2020); they find that porous aggregates have a higher chance of growing than compact aggregates. Thus, they can decouple from the gas present in the disk, which lets them survive and allows the possibility of further growth to planetesimal sizes.

Binary collisions between dust aggregates have been studied over the last 15 yr using granular-mechanics codes. Wada and coworkers have studied the influence of aggregate porosity, the collision velocity, and the impact parameter on the collision outcome of equal-mass collisions (Wada et al. 2007, 2008, 2009). They also investigated sequential collisions of aggregates and the induced aggregate compression (Suyama et al. 2008, 2012) as well as the outcome of mass-asymmetric collisions (Hasegawa et al. 2021). A key finding of these latter studies is that aggregate porosity is not a fixed characteristic of aggregates but evolves during their collision history. Using the same code, Seizinger & Kley (2013) and Seizinger et al. (2013) focus on aggregate growth and the bouncing threshold. This group also used smoothed-particle hydrodynamics simulations of aggregate collisions (Geretshauser et al. 2010; Meru et al. 2013), which allow the simulations to be extended to larger aggregates (on centimeter and decimeter scales); they report that aggregate porosity is key to the ability of aggregates to survive collisions and grow.

In this work we first present in Sect. 2 a brief summary of our previous methods and results that will serve as support for this study. In Sect. 3, we describe in detail the processes carried out by our Monte Carlo code, which we call the Collision of Porous Aggregates (CPA) code, starting with the construction of the aggregate population until the final mass, size, and porosity distributions are obtained, after a specified number of binary aggregate collisions occur. In Sect. 4, we apply our code to several representative cases and present the obtained results. We focus on a relatively small number of consecutive collisions that represent the early stages of aggregate growth, well below the gravity-dominated limit, and are crucial for protoplanetary disk evolution. Finally, in Sect. 5, we provide some conclusions and discuss several future possibilities for the code.

2 Collision outcomes from granular mechanical simulations

2.1 Granular-mechanics algorithm

In the recent past, we carried out numerous granular-mechanics simulations between porous silica aggregates using the well-documented granular mechanics package LAMMPS (Plimpton 1995). All our aggregates were built based on the model published by Ringl & Urbassek (2012), where granular interactions are based on the work by Dominik & Tielens (1997). It assumes a repulsive normal force based on the Hertzian model, additionally viscoelastic dissipation of the motion in normal direction, an attractive force, sliding friction, rolling friction, and friction of twisting motion. Of particular importance is the attractive force, which we modeled according to the Derjaguin–Muller–Toporov model (Derjaguin et al. 1975; Maugis 2000).

The aggregates consist of a collection of spherical grains that have the same properties, including radius of grain (R_graįn = 0.76 µm), mass (m_graįn = 3.68 × 10⁻¹⁵ kg), and density (2 × 10³ kg m⁻³). These grains only interact with each other if the distance between their centers, d, is d < 2R_graįn. The length δ = 2R_graįn − d is called overlap, and the grains will interact only if δ > 0. A cluster is defined as a set of connected grains, where the distance between the centers of two connected grains must be ≤2R_graįn = 1.52 µm. These aggregates can be built with a specific filling factor (ϕ = 1 − porosity), shape and number of particles n (for more details, see Ringl et al. 2012b; Planes et al. 2021). Next, we summarize some previous studies that are important for the development of this work, all of them have used the same model mentioned above.

2.2 Background: Previous work

In Planes et al. (2017, hereafter Paper I), we explored the impacts of spherical granular aggregates (projectiles) composed of a variable number of grains 5 ≤ n_P ≤ 500 on a large cubic granular target composed of n_t = 70000 grains. Both the projectile and the target have the same value of filling factor (ϕ = 0.36). They were constructed using the method of Ringl et al. (2012b) by filling grains homogeneously into a box until the required filling factor is reached. We denote with µ to the mass ratio between the target and the projectile. As all grains have the same mass we can take µ = n_t/n_p, so in this work the range of mass ratio studied was 140 ≤ µ ≤ 1400. Initially, the projectile is set at a position above the target so that there is no interaction with it. Then the simulation is started by giving each grain in the projectile the same velocity, υ, which we vary between 5 and 200 m s⁻¹. We investigate the dependence of crater formation and grain ejection during the collisional processes on projectile initial velocity υ and size n_p (or equivalently, on the mass ratio µ).

In Planes et al. (2020, hereafter Paper II) we explored mass-asymmetric collisions between spherical aggregates with µ = 60 at υ = 100 m s⁻¹, with particular emphasis on their dependence on aggregate porosity. Even with this large value of mass ratio, the results were very different respect to impacts on a flat bed or compared with Paper I. For high filling factors (ϕ > 0.2), a crater is formed in the target. We find here that – in contrast to impacts on granular beds – the crater rim fractures, forming petals. For low filling factors, the projectile can simply penetrate through the target, leaving a large hole. The separation between the two regimes is quite abrupt and occurs for filling factors of ϕ ~ 0.2. In the window of 0.20 < ϕ < 0.35, we observe the target to grow by the assembly of mass from the projectile. At larger filling factors, grain ejection increases, leading to a net mass-loss, while at smaller ϕ, the projectile tears the target aggregate. We also analyzed the compaction of the granular material after the collision and found that, in general, the target remnant strongly compacted and the average number of contacts between grains has almost doubled from its original value. Our results also showed that in mass-asymmetric collisions, even at high velocities, aggregates may not completely fragment but may even grow.

In our last work Planes et al. (2021, hereafter Paper III), several impact velocities, porosities, and mass ratios between the aggregates were used to determine the threshold values that could separate the agglomeration from the fragmentation of the target aggregate. We can classify the result of simulations into three possible outcomes:

• Sticking by penetration (SP): The projectile hits and penetrates the target (Güttler et al. 2010).

• Total destruction (TD): This occurs when n_large/n_tot ≤ 0.5 (n_large is the number of grains of the largest fragment and n_tot = n_p + n_t) according to Wada et al. (2008).

• Two fragments (TF): This is a previously unobserved outcome. The collision results in two fragments: one of them has the structure of a hollow cylinder, and the other contains a part of the target remnant and the majority of the projectile, forming a structure similar to a hemisphere. This process results from the piston effect mentioned in Paper II.

Figure 1 shows a snapshot with examples of (a) SP, (b) TF and (c) TD outcomes (For more details, see Paper II and Paper III).

In Paper III, erosion and accretion efficiencies were also studied, and the minimum impact energy required to fracture the sample was discussed. Compaction and the mass distribution of the fragments produced after the collision had been analyzed.

3 Monte Carlo code

The Monte Carlo method is a numerical method for solving mathematical problems through the simulation of random variables (Sobol 1994). It is a computational algorithm with a simple structure configured to follow a random path, which can be repeated N_* times, each of these N_* paths being independent of the rest, and then the N_* results are averaged. The error of the method is inversely proportional to $N_{*}^{1/2}$. Therefore, the more times the path is traveled independently, the more certain the result will be; however, there will always be a margin of error associated with this process. In our context, we started from an initial distribution of granular silica aggregates and used this method to simulate a random collisional chain between them. A collision chain is understood as the successive process of collisions in pairs that the population aggregates have. Then, a large enough number of independent chains are averaged to obtain the final distribution of resulting aggregates. The result of each individual collision that occurs between two aggregates will depend on the results obtained in our previous works (see Sect. 2), where the size, mass and porosity of the resulting aggregates were obtained according to the initial conditions of the collision. A summary of the variables used in the code and the text is included in Table 1.

Fig. 1 Snapshots at 50 µs showing the collision outcome – SP (sticking by penetration), TF (two fragments), or TD (total destruction) – between aggregates with µ = 10 and: (a) ϕ = 0.15, υ = 10 m s−1, (b) ϕ = 0.15, υ = 20 m s−1, and (c) ϕ = 0.40, υ = 50 m s−1.

3.1 Construction of the population

3.1.1 Input parameters

The code will receive as input a set of initial parameters, which can be divided into three groups. The first is the parameters related to the initial population of aggregates:

• Initial density (N_agg,0): This parameter indicates how many aggregates are generated at the start of the simulation.

• Initial size distribution: A distribution of initial radii R, (F(R)₀), is required as input. This distribution includes N_agg,0 initial aggregates, having R between a minimum value R_low and a maximum value R_high.

• Initial filling factor distribution: As for radii, an initial distribution of filling factors ϕ, is required, as (F(ϕ)₀). The N_agg,0 initial aggregates have ϕ values between a minimum value ϕ_1οw and a maximum value ϕ_high.

The second group comprises parameters related to the individual collisions that will occur:

• Impact velocity (υ): Initial velocity of all the individual grains that make up the projectile (smallest aggregate of the pair of aggregates to collide).

• Impact parameter (b): the impact parameter is the distance between the direction of the projectile’s velocity and the center of mass of the target.

Finally, the third group comprises parameters related to the collisional process:

• Number of collisions (n_co11): This parameter indicates how many collisions will occur consecutively in a collisional chain.

• Convergence cutoff (C_conv): This parameter indicates the expected statistical error, when this value is reached the simulation ends. It will be explained in detail later.

• Maximum number of chains (N_maxchain): This indicates how many chains will be executed at most, if C_conv is not reached the code will run until the N_maxchain chain.

For this work, we started all our simulations with N_agg,0 = 10⁵ aggregates. For F(R)₀ and F(ϕ)₀ we have distributions with features shown in Table 2.

The collisional evolution of the dust population requires some velocity distribution that is sampled by the CPA code. This velocity distribution in general is not well known and might even change with time. For example, in some astrophysical environments, dust-carrying gas will be turbulent, including temporally and spatially varying accelerations from eddies. Some previous works have dealt with determining relative collision velocities in these cases (Ormel & Cuzzi 2007). For certain scenarios, a Maxwell-Boltzmann (MB) distribution has been used (Windmark et al. 2012; Drazkowska et al. 2014), while for several cases it has been pointed out that a distribution with a fat tail, compared to an MB distribution, would be required (Windmark et al. 2012, and references therein), for instance a Levy distribution (Garaud et al. 2013). Velocity distributions are often associated with temperature. However, it is important to point out that temperature is sometimes not a well-defined quantity for granular materials (Brilliantov et al. 2018; Baule et al. 2018; Puglisi et al. 2017). In addition, a kinetic temperature associated with a velocity distribution of the grains will not be related to the temperature associated with radiation emission properties since these are determined from atomic rather than granular properties. The role of the chosen velocity distribution has to be carefully evaluated for each astrophysical scenario where the code would be applied. In this paper, for the sake of simplicity, we sample granular aggregate collision velocities from an MB distribution with υ_rms = 10 m s⁻¹. Different υ_rms can be readily obtained, and different distributions could be relatively easily added to the code. For instance, the current version also includes a gamma distribution, which might be relevant for scenarios where the high velocity tail would be thinner than for an MB distribution.

The influence of the impact parameter b on granular collisions has not been analyzed by us, since all the collisions we carried out have been central. Therefore, the influence of this parameter will not be taken into account in this current work.

Inputs of group (3) will be variable. We discuss the convergence of the code according to the results for different n_coll later.

3.1.2 Building initial aggregates

In this section, the initial N_agg,0 aggregates are constructed as follows:

1. A size, R_i, is randomly chosen from the distribution F(R)₀.

2. A filling factor value, ϕ_i is randomly chosen from the distribution F(ϕ)₀.

3. The number of grains in the aggregate, i (n_i), is calculated, assuming a spherical shape: $$n_i=\frac{4\pi {\varphi }_i{R}_i^3}{3V_{\text{grain}}},$$(1)

   where V_grain is the volume of one grain (V_grain = 1.8388 × 10⁻¹⁸ m³ for silica grains).

4. This aggregate is saved and the construction proceeds in the same way until N_agg,0 aggregates are obtained.

5. The total number of individual grains in the N_agg,0 aggregates, nagg,tot, is calculated.

Then the initial distribution of mass depends on the initial size and porosity distributions. In this work, distributions of Table 2 give an initial mass distribution that exhibits power-law relationship with exponent ≃−2/3.

3.2 Selection of aggregates to collide

In this section, we will detail how the 2 aggregates that will collide with each other are selected from the total aggregate population. It is important to note that as the n_coll collisions occur, the total number of aggregates (N_agg,j) can vary after each jth collision (where j = 1, 2,…, n_coll), but by conservation of mass, the total number of individual particles in all aggregates nagg,tot, must be constant at all times. The values in the size (filling factor) distribution start at R_low (ϕ_low) and increase by steps whose width is R_step (ϕ_low), with values given in Table 2. There is no a priori maximum for the size distribution, but for the case of the filling factor distribution the maximum final value is 1. The code rounds the numbers obtained after the collision to the appropriate binning value. Using a finer binning than the one indicated in Table 2 did not affect our results, but with a coarser binning some features were lost.

The collision probability between two aggregates of radius R₁ and R₂ is related to their cross-sectional area. Therefore, the collision “kernel” must be proportional to (R₁ + R₂)². Other kernels have been used in the past, and analytical models are available for simple kernels, when only coagulation is considered. Results for a constant kernel where aggregates are chosen at random, independently of their size and velocity, and for a linear kernel where the choice is related to the sum of aggregate masses, are shown in Appendix B.

Here again, we emphasize that no study has been conducted on the result of collisions where the aggregates to collide have different porosities. Therefore, we do not have information to incorporate this type of collision into our code, and this type of collision is not allowed in this work.

In this work the selection of aggregates follows the direct simulation Monte Carlo technique developed by Bird (1963):

1. For collision j in a given chain, two aggregates, A₁ and A₂, are randomly selected with the same filling factor from the size distribution.

2. A velocity, υ, is randomly chosen from the given velocity distribution, F(υ). This is be the impact velocity between A₁ and A₂. F(υ) has a maximum velocity υ_max.

3. The code calculates a collision probability using a collision kernel proportional to the cross section of A₁ and A₂. $P_{\text{coll}}=\left[\upsilon \cdot {\left(R_1+R_2\right)}^2\right]/\left[{\upsilon }_{\max }\cdot \left(R_1+R_{\max }^2\right)\right]$, where (R₁+ R₂)_max is the sum of the two largest sizes in the aggregates population.

4. A random r value is taken from a uniform distribution between 0 and 1.

5. If P_coll > r, the collision between aggregates is performed; otherwise, the code repeats the process and selects two new aggregates, starting again from the first step.

6. The mass ratio, µ, between A₁ and A₂ is calculated. Since all individual grains have the same mass, µ can be calculated as the ratio between the number of individual grains, n_i, that make up each aggregate A_i.

7. According to µ and υ, the result of the collision is determined (see the next section). The resulting aggregates are incorporated into the distribution, which is saved as N_agg,j+1.

8. After the collision is performed, the two aggregates A1 and A₂ are removed from distribution N_agg,j.

9. It is verified that n_agg,tot is the same for both N_agg,j· and

10. The process is repeated from the first step, with j = j + 1, until j = n_coll.

Therefore, each individual collision will have the following parameters:

• Aggregate A₁: R₁, ϕ, n₁

• Aggregate A₂: R₂, ϕ, n₂

• Mass ratio, µ

• Impact velocity, υ.

3.3 Possible results after a single collision between two aggregates

The CPA code randomly chooses two aggregates and their impact velocity, only imposing that they have the same filling factor. In this section we assign a result to each of the individual collisions. First, we classify its outcome as SP, TF, or TD – definitions according to Sect. 2 – and then we obtain the result of the collision: it must include the number of resulting fragments, the number of particles in each of them and their porosity.

3.3.1 Binary collision outcomes

Numerically, it is not feasible to run all the possible combinations of µ, ϕ, and υ to obtain the exact information about the aggregates resulting in each collision. Therefore, based on our previous results, summarized in Sect. 2, we divided our space of parameters according to cutoff values.

First branch, mass ratio: Some previous works found different results when both collision partners had the same or different mass (µ = 1 or µ ≠ 1). In Paper III we find that a single cutoff value seems not to be enough, so as a first improvement (and approximation) we divide the results according to four possible intervals of µ (the exact values that we previously analyzed are: µ = 1, 10, 60, and 140 < µ).

Second branch, filling factor: The outcomes for the different collisions have been separated between very porous aggregates (ϕ < 0.20) and porous aggregates (ϕ > 0.20). In some cases, it was necessary to add an intermediate interval that exhibited a different behavior in our previous simulations (0.20 ≤ ϕ < 0.30).

Third branch, impact velocity: The last, but no less important, parameter to take into account to separate the regimes is the impact velocity, υ. In some cases we have precisely defined the critic impact velocity, υ_crit, that separates the outcome SP from the TF or TD (Paper III). But in some cases we have only a velocity range where the transition between regimes occurs: in these cases we made an estimate of υ_crit.

Based on this, we established ranges in each parameter to classify the collision as SP, TF, or TD, as shown in Table 3 with their corresponding references. A visual outline of the possible results for each collision can be found in Fig. A.1.

3.3.2 Resulting fragments

The code has so far classified the result as SP, TF, or TD. We describe below what implications the different regimes have. We remember that from each collision we want to obtain the total number of fragments that result and, for each of these resulting fragments: the number of grains that compose it and its final porosity (also its size, but its radius can simply be calculated with these last two values, taking Eq. (1) into account).

After the collision, one or some large fragments can remain and the rest of the particles will be part of the ejecta (very small fragments). For the size distribution of the resulting fragments we have the following parameters:

• n_tot: The total number of particles in both aggregates, projectile and target (n_tot = n_p + n_t).

• Y_s: The number of particles ejected in small groups or dust. Based on our previous DEM simulations (Papers I, II, III), the maximum size for the resulting fragments was set to 100 grains for this work.

• Y_L. The number of particles ejected into large fragments (with more than 100 grains).

For example, Y_s does not include the grains in the largest aggregate result in SP, or in the two largest aggregates resulting in TF. For TD the situation is different: when υ is high enough, all the ejecta is dust (Y_s ≃ n_tot), but when v is close to υ_crit, we observe a number of large fragments (between 1 and 30), and the rest is dust. This detail will be taken into account later.

Based on our results from Paper III and references here, we use a power-law relationship to represent the size distribution of the particles ejected in small fragments: $$F\left(n\right)=k{n}^{-\tau },$$(2)

where k is a constant. So we need to know the value Y_s: how many of the total grains will be part of the ejecta after the collision, and the value of the exponent τ. In Table 3 we present the results for this value obtained in previous works (when not detailed, it means that these data were not measured). Taking this into account, we assigned characteristics to each outcome in order to maintain a simple model that can be efficiently executed by our code, but that at the same time it satisfactorily represents all the results previously obtained by us. As always, any of these values can be changed by the user if other criteria are required. For each outcome our code adopts the following Y_s, Y_L, and τ values. For SP, Y_L = 0.99n_tot (a single large fragment) and Y_s = 0.01n_tot with τ = 2.8. The only exception is µ > 100 where Y_s = 100 and all of them are monomers. For TF, if 1 < µ ≤ 20: = 0.98 n_tot (Fragment 1: 0.49 n_tot; Fragment 2: 0.49 n_tot) and Ys = 0.02 n_tot with τ = 2.8. If 20 < µ ≤ 100: = 0.98n_tot (Fragment 1: 0.84n_tot; Fragment 2: 0.14n_tot) and Y_s = 0.02n_tot with τ = 2.8. We remark that Y_L has a strong dependence on υ (see Table 3 and references there). Therefore, we took only two possibilities to simplify the model. Finally, for TD, Y_L = 0.7n_tot (five equal-sized fragments) and Y_s = 0.3n_tot with τ = 2.8. If there are initially n_tot ≤ 5 particles, after the collision they will be all monomers.

In the TD case, Y_L and Y_s also depend strongly on υ. For example, according to Table 3, Y_s varies between 0.2n_tot (for υ close to υ_crit) and 0.99n_tot (for the highest velocities, υ, analyzed in this work).

At this point, two clarifications are necessary: (a) in both TF and TD cases, it has been preferred to take the values of Y_s and Y_L that are observed when the impact velocity υ is close to υ_crit. This is because the v distribution is an MB distribution with mean of 10 m s⁻¹. Therefore, if a simulation results in TD, for example the case µ = 10, ϕ = 0.15 in Fig. A.1, it will be more likely to occur at υ = 35 m s⁻¹ than at υ ≪ 35ms⁻¹. (b) The code rounds the number of particles in each resulting fragment to integers. If, for this reason, the final number of particles after the collision exceeds the value n_tot, the surplus grains will be subtracted from the monomers. If this number is less than n_tot, the missing grains will be added as monomers.

3.3.3 Porosities of the resulting fragments

As aggregates are very porous, during collisional process they could restructure and change their filling factor. The compaction factor is defined as the ratio of the filling factor ϕ before and after the collision.

We have limited data on compaction changes after a collision, but we found that it depends on µ (Paper III): if we consider a υ = 10 m s⁻¹ (υ_rms value of our MB distributions of random initial velocities), for all the φ analyzed, when µ = 10 we observe an increase in compaction of 75–100% and when µ = 60 of 50–75%. Gunkelmann et al. (2016b) find for µ = 1 (and 0.08 < ϕ < 0.201) an increase of 100%. So as µ increases the final compaction decreases. From Fig. 5 in Paper III, and Fig. 7 in Gunkelmann et al. (2016b) we see that a 100% increase seems to be approximately the maximum reached and then the compaction values decrease.

Because of that, we assigned two different compaction ranges: between 1.5 and 1.75 if µ > 20 and between 1.75 and 2 if µ ≤ 20. According the value of µ, the code will randomly choose a compaction factor from the corresponding range. Of course, this is a coarse approximation, since compaction depends on other factors, for example, from Fig. 5a in Ringl et al. (2012a) we note that the central collisions (like our simulations) present a maximum in the compaction factor. But in the absence of more information in this regard, this is an improvement over considering a single compaction factor for all collisions, and can still be easily modified in the code when there are more simulations to support better regime discrimination.

We set a maximum limit for this value at ϕ_max = 0.74, which is the maximum possible packing for spheres of the same size $(\pi /\sqrt{18}=0.7404$, considering a face-centered cubic structure; Bargiel & Tory 1993). Although some authors suggest that this value is not reached by collisions that cause disorderly compaction, suggesting maximum values of ϕ_max ~ 0.58–0.65 (Torquato et al. 2000; Weidling et al. 2009; Teiser et al. 2011). We emphasize that our objective is to evaluate the global compaction of these aggregates, then, we prefer to cover the largest filling factor theoretically possible, which can then be limited in the post-analysis.

There are two important exceptions: The first are monomers: aggregates that have a single particle (n = 1) will have an assigned filling factor of ϕ = 1. The second are dimers: aggregates made up of n = 2 particles should have an assigned filling factor of ϕ = 0.25 (two touching spheres radius R_graįn enclosed by a sphere of radius 2R_graįn), although the code assigns them a different value, ϕ = 0.9125, in order to distinguish them from aggregates that have ϕ = 0.25 and n ≠ 2. This value can then be filtered or changed in the post analysis.

This will allow a cutoff to be put at ϕ ~ 0.75 to observe the porosities of aggregates with n ≥ 2 since the porosity of aggregates with few particles is meaningless.

3.4 Final distributions

3.4.1 Resulting distributions after n_coll collisions

At this point, the code executes j = n_coll collisions, following the procedure described in Sect. 3.2 and using the criteria set out in the previous section for each individual collision. We want to get: the mass (F(n)), size (F(R)), and porosity (F(ϕ)) distributions of the final resulting aggregate population. Furthermore, obtaining a matrix (F(n, ϕ)) that interrelates the mass and porosity data of the aggregates can be interesting.

Then, once the last collision ends, the code returns each distribution normalized with the total number of aggregates existing in the final population $N_{\text{agg,}{n}_{\text{coll}}}=N$:

Mass distribution, F(n). This is the number of aggregates composed of n grains. Since the mass of all grains is identical, this can be considered a normalized mass distribution with m_graįn. Here, the possible values of n are hundreds of thousands (and they will rarely match on all chains to simulate later). Therefore, a grouping has been made in logarithmic bins.

Size distribution, F(R). This shows the radius of all the final aggregates grouping them into intervals that are 1 µm wide (the width is double the initial to smooth the curve). Each frequency is located in the initial value of this interval.

Porosity distribution, F(ϕ): It is presented through the distribution of filling factors, which groups the resulting values into intervals whose width is 0.05 (the width is double the initial to smooth the curve). Each frequency is located in the middle value of this interval. Remember that, by default, the code assigns a value of ϕ = 1 to monomers and ϕ = 0.9125 for dimers.

Mass-porosity distribution, F(n, ϕ): This matrix details the frequency of aggregates with a certain mass (grouped in logarithmic bins) and filling factor values, making it possible to observe the porosity of the resulting aggregates according to their mass. It also includes a column with the respective size for the n bin (where Eq. (1) is used) and so also estimates the size-porosity distribution (F(R, ϕ)). This last distribution considerably increases the computation time, and it is the user's choice to request it.

3.4.2 Resulting distributions after N_chain chains

So far we have detailed the results after n_coll collisions occurred in one chain (one path traveled). Now we repeat this process starting from the same initial population of N_agg,0 aggregates, but new values will be randomly chosen to repeat the process detailed in Sect. 3.2, N_chaįn times, where this number must be large enough to guarantee the convergence of our results. Next we detail how the total number of chains to be run, N_chaįn, is defined according to the desired error, and the calculation of the final distributions averaged over N_chaįn chains. We also estimate the standard deviation of each point for these distributions.

We distinguish two types of errors: (a) bin error is the calculation of the standard deviation of each bin of each distribution and (b) global error is the global error of the code. To explain this procedure, we need to define some new quantities: The convergence cutoff, C_conv, indicates the accepted global error; once reached, the code will end. Finally, N_maxchain is the maximum number of chains. The N_chain is the number of chains that the code has traversed until reaching the convergence, C_conv.

For each distribution i normalized by N, F(i), the procedure is the same:

1. The value of each bin of F(i) is saved for the first chain that has been traversed, named F(i)^avg,0.

2. The next chain is simulated, and the value of each bin of F(i) is saved.

3. The average value of each bin of F(i) is calculated using the values from both chains: F(i)^avg,1.

4. The equation Δ(i) = F(i)^avg,1 − F(i)^avg,0 is calculated.

5. Only one of the distributions, i, is used to calculate the convergence of the simulation (i.e., i = n, R or ϕ). For the selected distribution, the global error is the sum of Δ(i)² over all bins.

6. If the global error is greater than the C_conv cutoff, F(i)^avg,1 is saved, F(i)^avg,0 is deleted, and the next chain is simulated, since adding new chains will reduce the error in the limit of large numbers of chains.

7. If the error is less than the C_conv cutoff during a few N_cut consecutive chains, each of the final distributions is calculated until the last chain (N_chain), $F{\left(i\right)}^{\text{avg,}{N}_{\text{chain}}}$, then the data are output to the disk, and the simulation is ended. We set N_cut = INT(1000/n_coll) for this work.

Any of the distributions can be chosen to calculate this. In particular, we chose F(R) since it obeys being the least grouped distribution and its data have not been modified (as is the case with some filling factor values). And once the convergence process is reached it is normalized by the total number of chains traveled, N_chaįn. Finally, each distribution is normalized to area 1, that is, Σ_i[F(i) · bw_i] = 1, where bw_i is the bin width; therefore, this normalization will depend on the width of the bins adopted in each case.

This code was tested by varying the seed, that is, generating a new population of N_agg,0 aggregates for each chain, and repeating the entire collisional process. No quantitative changes were observed either in the initial population generated (of 10000 aggregates) or in the final distributions. Therefore, we conclude that it is not necessary to modify the seed in each chain of collisions traveled. Again, this is left to the choice of the user.

To summarize, the final distributions originate not from a single fragmentation event, but from multiple ones, given by n_coll. To obtain statistically meaningful results, N_chain collision chains have to be averaged to obtain a “steady” state, where the simulation of an additional chain will not affect the average results beyond the limit given by the convergence parameter. Therefore, the distributions are “steady” with respect to N_chain (keeping n_coll constant), not with respect to a possible increase of n_coll. CPA users can modify this convergence parameter. As expected, more collision chains are required as the average number of collisions experienced by aggregates decreases. The CPA code is shown to converge reasonably well for the scenarios studied here, and can thus be run in standard desktop workstations.

The complete flow diagram of the Monte Carlo software is shown in Fig. 2.

Fig. 2 Flow chart diagram of the Monte Carlo simulation.

Fig. 3 Mass distribution (normalized by area). Colors denote the initial distribution and final distribution after 5, 10, and 100 collisions. Gray lines show fits for the exponent α (see text).

4 Results

In this section, we show the results of the CPA code. We ran the simulations with the input parameters detailed in Table 2, C_conv = 10⁻⁶ and three values for n_coll. 5, 10 and 100. We only present some representative results for each of the distributions. We focus on cases with relatively low n_coll, which are relevant for regions with very low dust density, and for the early stages of dust aggregation, well below the gravity-dominated regime.

4.1 Mass distributions

Figure 3 shows mass distribution, F(n), of the initial population (black color) and of the final population after 5, 10, and 100 random collisions. This histogram represents the frequencies of the number of grains in each aggregate, n_i, grouped into bins [i, i + w), where the widths of intervals w are logarithmic. The initial distribution of mass (created following Sect. 3.1.2) presents a power-law relationship: $$F\left(n\right)=a{n}^{\alpha },$$(3)

where α is a constant and α = −2/3. After n_coll collision this value for α holds for aggregates with n > 30 when n_coll = 5, 10 but only for aggregates with n > 100 when n_coll = 100, as seen in the fit in Fig. 3. For the smallest aggregates a different exponent α is observed: −2.8. This is reasonable, because this is the τ value that was assigned to the resulting small fragment size distribution. As more collisions occur (higher n_coll), a greater number of fragments will exhibit this behavior; for this reason, when n_coll = 100, the exponent τ holds up to a greater n (100). Preliminary results for the full model with much larger n_coll show a steeper distribution, with larger α, as expected from the models by Birnstiel et al. (2012a). Some results for simple coagulation kernels and n_coll = 10000 are shown in Appendix B, showing excellent agreement with analytical results.

4.2 Size distributions

From Table 2, our initial size distribution is uniform, as we can see in Fig. 4 in black color. As in the case of mass distribution, the flat distribution holds for large aggregates, and other power law is observed for small aggregates: $$F\left(R\right)=b{R}^{\beta },$$(4)

where b is a constant. A value of β = −6 is shown in Fig. 4. This fits works well for aggregates until R ~ 5 µm when n_coll = 5, 10 and until R ~ 10 µm when n_coll = 100. Again, this makes sense: as more collisions occur, more aggregates are generated and leave the initial uniform distribution. Another interesting thing is the generation of aggregates with R > R_max, it means, at the end, the population contains larger aggregates that in the beginning. These larger aggregates have a low frequency compared with aggregates with R < 100 µm. However, as n_coll increases, this difference in frequency decreases as well. Of course, Fig. 4 shows the global results, and really large aggregates can appear: our results include R until 216 µm, but these cases have a frequency that is too low and can be neglected. We note that this slope is quite different to the one used to fit ISM dust size observational results (Mathis et al. 1977). Collisions and porosity are not expected to play a dominant role in that scenario.

Sometimes it is possible to relate the power-law exponent for F(R) to that for F(n) found in the previous section. The number of aggregates dN with radii between R and (R + dR) is defined by the differential size distribution: dN = N(R)dR = C_RR^sdR, and with masses between m and (m + dm) by the differential mass distribution: dN = N(m)dm = C_mm^qdm, where C_R and C_m are constants. Then it is $$\text{d}N=C_R{R}^s\text{d}R=C_m{m}^q\text{d}m.$$(5)

For spherical aggregates, it is $$m=\frac{4}{3}\pi \rho \varphi R^3.$$(6)

Using Eq. (6) in Eq. (5), we obtain $$C_R{R}^s\text{d}R=C_R{\left(\frac{3m}{4\pi \rho \varphi }\right)}^{s/3}\text{d}R.$$(7)

Assuming ρ and ϕ constants and differentiating Eq. (7) yields $$\text{d}m=\frac{4}{3}\pi \rho \varphi \left(3R^2\right)\text{d}R\Rightarrow \text{d}R=\left(\frac{\text{d}m}{4\pi \rho \varphi }\right){\left(\frac{4\pi \rho \varphi }{3m}\right)}^{-2/3}.$$(8)

Replacing Eq. (8) in Eq. (7), and after some algebra, one obtains $$R^s\text{d}R=\left(\frac{1}{4\pi \rho \varphi }\right){\left(\frac{3}{4\pi \rho \varphi }\right)}^{\frac{s-2}{3}}{m}^{\frac{s-2}{3}}\text{d}m.$$(9)

This allows one to obtain the relationship between the exponents s and q as $$q=\frac{s-2}{3}.$$(10)

In our simulations, the filling factor is not constant, but we can still test these approximations. We use α and β for the exponents in the numerical fit to our results for mass and size distributions, respectively. For our initial size distribution s = 0 and we obtain α = −0.66 ~ q = −2/3, as expected from Eq. (10). For the small aggregates obtained after the collision chains, we fit a power law for the size distribution, with β = −6, which would be equivalent to use s = −6 in Eq. (7), giving q ~ −2.67 from Eq. (10). We obtain a value of α = −2.8, fairly similar to the analytical approximation. The main reason for this agreement with the analytical approximation is the relatively low number of collisions studied here, such that the mean value of the filling factor would remain roughly constant.

In our simulations, the initial filling factor of each aggregate is taken randomly from a uniform distribution, but the code does not provide the evolution of the distribution of average filling factor versus aggregate size. The relatively low number of collisions explored in this work is expected to cause only slight changes to that distribution. Future work exploring longer collision chains might explore this further, as it has been done before in coagulation simulations (Ormel et al. 2007; Okuzumi et al. 2012).

Fig. 4 Size distribution (normalized by area). Colors denote the initial distribution and the final distribution after 5, 10, and 100 collisions. Gray lines show fits for the exponent β (see text).

Fig. 5 Porosity distribution (normalized by area). Colors denote the initial distribution and the final distribution after 5, 10, and 100 collisions. (a) Complete distributions. (b) Only aggregates with n > 2.

4.3 Porosity distributions

Figure 5 shows the initial and final filling factor distributions for (a) all aggregates and (b) aggregates with n > 2. Figure 5 (a) is in agreement with our previous plots: monomers (with a filling factor fixed in ϕ = 1) are predominant, and dimers are in second place, respecting the pre-accorded filling factor of ϕ = 0.925 (values according to Sect. 3.3.3). But talking about filling factor of monomers or dimers has not physical sense. So in Fig. 5 (b), we remove these aggregates to focus in aggregates with n > 2 (it is a zoom of (a)). Remember that the value ϕ = 0.75 was chosen as the upper limit for the compaction, so this figure shows that as n_coll grows, more aggregates find this upper limit as a consequence of collisional processes. Also, for n_coll = 5, 10 the rest of the distribution seems to hold quite uniform, but for n_coll = 100, the distribution seems to take an increasing trend in ϕ.

4.4 Mass and size distribution for different ranges of aggregate porosities

For this analysis, the possible porosity values will be divided into three groups: (1) P (very porous): ϕ ≤ 0.2, (2) PC (porous): 0.2 < ϕ ≤ 0.4 and (3) C (compact): ϕ ≥ 0.4. The separation between PC and C follows the one proposed by Güttler et al. (2010), and between P and PC is due to our detailed results in Paper II. For each interval in ϕ the corresponding mass and size distribution will be presented. This proposed representation is our choice, the code gives freedom so that other alternatives can be proposed according to the object of its application. The purpose here is understood if the previous distribution shown for all aggregates present any dependence according the aggregate porosity values.

4.4.1 Mass distribution according to porosity ranges

Figure 6 shows the mass distribution, in logarithmic bins, for the 3 intervals of porosity for: n_coll: (a) 5, (b) 10, and (c) 100. We indicate Eq. (3) with α = −0.66 and −2.8 as in Fig. 3. Both exponents fit quite well. However, the group C covers the largest interval of values in n. We note that monomers and dimers only belong here, so C starts at n = 1 whereas PC and P at n = 3. On the other hand, the largest aggregates are produced by a collisional process, where compaction took place, so PC covers higher values of n than P, and C than PC. Beyond the α’s values being maintained, frequency in smaller aggregates increases as n_coll increases. Erosion and generation of small clusters are an extremely important part of the collisional process.

Appendix C shows distributions separated by porosity region.

4.4.2 Size distributions according to porosity ranges

Figure 7 shows the size distribution for the 3 intervals of porosity for: n_coll: (a) 5, (b) 10, and (c) 100. This information for R was obtained throw Eq. (1) (see Sect. 3.4.1), but F(n) is grouped in logarithmic bins, so an error associated with this grouping must take into account. We show these results as an example, the user can request that the code calculate F(R, ϕ) directly, but F(R) has ~220 bins and F(n) has ~60, so it will be more expensive in computational time.

Using Eq. (1), Fig. 7 shows dependences dragged from F(n). We plot two values for the β exponent of Eq. (4): −6 as in Fig. 4 and −2.8 as in Fig. 3. By the way the aggregates are built, now the three regimes, P, PC and C start at R = 1 µm, but again, larger aggregates are more compact. We observe that for small aggregates the relationship proposed (with β = −6) fits very well the cases n_coll = 5 and 10, but for n_coll = 100 the slope seems to be steeper. Class C includes mostly new aggregates generated by collisional processes. The distribution of F(R) has no “grouping” into logarithmic bins to smooth the curves, as F(n) does. For this reason and for a relatively low number of collisions where new values of R appear with low frequency, there is a large scatter in Fig. 7.

Fig. 6 Mass distribution (normalized by area) for ncoll: (a) 5, (b) 10, and (c) 100. Colors denote ranges of the porosity of the aggregates. Gray lines show fits for the exponent α (see text).

Fig. 7 Size distribution (normalized by area) for ncoll: (a) 5, (b) 10, and (c) 100. Colors denote ranges of the porosity of the aggregates. Gray lines show fits for the exponent β (see text).

4.5 Convergence at large radii

Convergence at large radii, larger than the maximum radius of the initial distribution, requires decreasing the convergence parameter C_conv· This implies a significant increase in computational time, and it might not be needed for some applications. As an example, we show changes for the case n_coll = 100 in Fig. 8 for the complete range of R (a) and for large radii (b), but similar behavior is observed for n_coll = 5, 10. For very larger R, as C_conv decreases, more continuity is observed, but remember that the low frequency makes these aggregates negligible. We zoom into the range 100 µm < R < 150 µm (Fig. 8 (b)), where the effect of decreasing C_conv is clear: dispersion is significantly reduced and for C_conv = 10⁻⁶ or lower there are no significant changes.

Fig. 8 Size distribution (normalized by area) for ncoll = 100 at different Cconv. (a) Complete range of R. (b) 100 µm < R < 1 50 µm.

4.6 Computational performance

Table 4 shows the total number of chains and the computational cost in the form of wall-clock time for a workstation with an AMD Ryzen 7 3800X 8-Core and 64GB DDR4 RAM at 2600MHz for C_conv = 10⁻⁶ and 10⁻⁸ as examples. Our code does not use a representative particle approximation, but can still model growth by binary collisions over many orders of magnitude thanks to significant hardware improvements since early pioneering MC work (Ormel et al. 2007), as also seen in Appendix B. Computational performance of the code is reasonably good, and up-scaling the runs to include a larger number of initial aggregates is possible, increasing the simulated range of aggregate growth. The maximum consumption of RAM memory for any of the cases shown in Table 4 is close to 4MB. Memory usage grows linearly with the initial aggregate number, N_agg,0. As an example of computational cost, a simulation with N_agg,0 = 10⁶, n_coll = 10⁶ and C_conv = 10⁻³ takes almost 26 hours and 150MB RAM for the constant kernel from Appendix B, in the same desktop computer used for the rest of the simulations. We also note that total running times increase linearly with N_agg,0. Assuming n_coll ~ N_agg,0 = 10^p one would obtain dynamic aggregate growth over p orders of magnitude for the kernels without fragmentation, shown in Appendix B.

Another issue is that the final number of aggregates for each simulated collision chain could be large, affecting computational performance. For example, the final number of aggregates averaged over all independent chains, N_agg, for an initial population with N_agg,0 = 10000, leads to N_agg = 10145 for n_coll = 3 and N_agg = 114 179 for n_coll = 1000.

Our code runs serially, in a single core. Future versions might add parallelization by sending different collision chains to different parallel threads, with a master process collecting the final distributions and normalizing them. This would ensure good parallel scaling if the number of chains, required for convergence is lower than the number of parallel processes. For the tests presented here N_chain ~ 10⁴−10⁷, and the future parallel code could perform in a supercomputer.

5 Conclusions and future possibilities

The CPA code was developed to study the collision processes between mass-asymmetric, spherical, and micron-sized porous silica aggregates that belong to a population of dust in a vacuum, without gravity. Our simulations, which use a relatively low number of collisions, n_coll, gave the following results.

1. The size distribution, F(n), is directly related to the mass distribution. For a distribution that is uniform in size, R, it is F(n) ∝ n^−2/3. As the population evolves due to collisions, this dependence remains unchanged at large sizes but is broken at small sizes, below a critical size, n_crįt, which increases as n_coll increases.

2. Collisional evolution produces a huge number of monomers and dimers, with small clusters following a power law in radius with exponent ~−6 (and ~−2.8 for a power law in mass), even for very short collision chains. Although it is not the main objective of these simulations, we show the feasibility of collisional growth of very large clusters. Studying further growth would require the inclusion of gravitational forces and other effects.

3. When collisional processes are assumed to end, larger aggregates could appear in the final distributions, in addition to a significant production of small aggregates (R < 10 µm). The frequency of the largest aggregates (R > 150 µm) is ~10 orders of magnitude lower than the small-sized aggregates, so they can often be disregarded.

4. Although porosity changes due to collisions are usually neglected, collisions lead to the compaction of the dust population, which is expected. As n_coll increases, a larger number of aggregates reach their lowest possible porosity because they are strongly compacted by one or more consecutive collisions. These results verify that a dust population subjected to a (relatively low) number of random collisions will present an important variation in terms of porosity distribution, which modifies the results of future collisions of the resulting aggregates (Paper II, Gunkelmann et al. 2016a).

5. In general, the porosity distribution, F(ϕ), shows a population at low porosity that increases with n_coll, along with the population for monomers (filling factor ϕ = 1) and dimers, which are given an artificially high ϕ (0.9) to separate them from the rest of the dust population. In general, larger aggregates result from agglomeration, while monomers and dimers arise from collisions involving fractures, with both processes usually leading to compaction.

6. The CPA code provides a matrix that contains the mass and porosity of all resulting aggregates. A relatively simple post-analysis can be carried out to separate mass and porosity into bins chosen by the users, to be examined separately. When F(n) is evaluated according to different ranges of porosity, the power law proposed for the entire distribution is roughly maintained. However, there are some differences for n_coll = 5, the smallest number of collisions considered here, which represents the very early stages of aggregate growth.

7. The CPA code calculates the size, mass, and porosity distributions for coagulation-fragmentation equilibrium after a certain number of collisions, n_coll. Statistical convergence of the distributions can be explored by users varying the “Convergence” parameter.

8. In the present study, we focus on the final distributions after a relatively low number of n_coll. However, one could imagine a time-dependent scenario where the number of collisions increases with time and the collision time between aggregates is given by their cross section and the velocity distribution. Therefore, one could equate n_coll with the elapsed time and study time evolution. This would require recalculating average collision times for each new size distribution.

This code incorporates several improvements over the features in similar previously developed codes (Birnstiel et al. 2012b; Zsom & Dullemond 2008; Ormel et al. 2007). From a methodological point of view, the stochastic collision model is based on the well-known and accepted direct simulation Monte Carlo approach, and an internal convergence check is implemented. From a physics viewpoint, the code emphasizes the inclusion of mass-asymmetric collisions (covering a wide range in aggregate radius), inter-granular friction, and the influence of porosity. Hereby, the code design follows a multi-scale approach that allows handshaking with results from granular-mechanics codes to provide input parameters on aggregate collisions. However, there is still much to improve. We emphasize that any improvement would likely entail an increase in computational time, both in the simulation time of the code itself and in studies or simulations that support new results to be incorporated. We believe that with current technological advances, the improvements proposed below could be achieved in the near future:

• Aggregates of different materials: In addition to silica, there is a need for simulations that provide information on collisions between aggregates of other materials, such as ice – of particular interest in astrophysics granular collisions (Wada et al. 2009; Kimura et al. 2020) – and organic (carbon-based) matter (Nietiadi et al. 2020b). Aggregates whose composition is a mixture of different materials (such as a silica core and a water shell) would also need to be considered (Nietiadi et al. 2020a).

• Porosity: There are two strong points to improve here. Firstly, simulations between aggregates with different porosities have not been studied extensively, and therefore this is not incorporated into our code. Nevertheless, results from a binary collision may change if the aggregates have different porosities, especially in extreme cases, for example, if a very compact aggregate impacts a very porous one. This would probably emphasize the piston effect observed in Paper II. Alternatively, if a very fluffy aggregate impacts a compact one, the projectile could perhaps be disarmed almost without disturbing the target. Despite these assumptions, further molecular dynamic simulations are required to test these hypotheses, and we plan to perform them in the near future. Second, there is no consensus on how the final porosity varies as a function of impact velocity and mass ratio. We have generally studied compaction as a function of the coordination number, which is related to the filling factor, φ. Nonetheless, studies focusing on the evolution of φ are necessary.

• Impact parameter: Although the impact parameter, b, is included as a variable in our code, we do not have enough information about the influence of this factor, and, therefore, it is not considered. However, based on some previous studies (Ringl et al. 2012a; Wada et al. 2013), growth is favored over erosion for central collisions. It would be interesting to be able to include the dependence of these results on b.

• Impact velocity distribution: The influence of the velocity distribution needs to be investigated. The choice of this distribution depends on the particular application of the CPA code, but even in the same scenario there seems to be no consensus about which distribution should be used (Windmark et al. 2012; Drazkowska et al. 2014; Garaud et al. 2013). In principle, in a more self-consistent code, the velocity distribution may need to be modified as the system evolves over time. In future work, we will analyze the changes in the final distributions of mass, size, and porosity when different velocity distributions are employed.

• Grouping of outcomes according to impact velocity: for simplicity, some regimes were separated according to a specific value in velocity. This may not be extremely precise because it was based on a limited number of granular simulations. New simulations covering more υ values for each set (ϕ, µ) are required in order to refine these boundaries. In addition, the number of resulting large fragments and their masses show a strong dependence on υ. However, we believe that it is possible to improve these points gradually as more results are obtained by different research groups.

• Aggregate form: All our aggregates are considered spherical, relating R and n. Although it is not easy to consider irregular geometric bodies, because that would require a huge number of collisions for different rotations of the aggregates, we are aware that this would be a much more realistic case for representing dust aggregates.

• Grain radius, R_grain: Results can vary when considering a different R_grain or a distribution of grain sizes. Umstätter & Urbassek (2020), for instance, show that bouncing was strongly enhanced for a bimodal grain size distribution. Compaction and fracture are also expected to change significantly (Pál et al. 2016).

• Number of collisions, n_coll: Future simulations could explore in detail the regime of much larger n_coll, allowing one to compare simulation results with the analytical results from Birnstiel et al. (2012a).

We expect to apply the CPA code to particular astronomical scenarios, but we have made the source code freely available in the hope that it can help other astronomy research projects that deal with porous dust collisions¹.

Acknowledgements

B.P., E.M., and E.M.B. acknowledge support from ANPCyT PICTO-UUMM-2019-00048 and SIIP 06/M008-T1. This work used the Toko Cluster from FCEN, UNCuyo, which is part of the SNCAD, MinCyT, Argentina.

Appendix A Summary of final binary collision outcomes

Fig. A.1 presents the classifications made by the CPA code for each binary collision, according to impact velocity, mass ratio and porosity of the aggregates, in accord with the information contained in the first four columns of Table 3. The definitions of regimens SP, TF and TD are found in Sect. 2.2, and how this classification affects the code is explained in Sect. 3.3.2.

Fig. A.1 Relation between initial conditions and final outcomes.

Appendix B Coagulation tests

Our CPA code is not designed exclusively for coagulation, but takes into account possible fragmentation and generation of monomers or small aggregates, together with porosity changes. Regardless, we tested our code with two common coagulation kernels (linear and constant), excluding fragmentation. Porosity and velocity distributions do not play any role for those two kernels. With the constant kernel, the aggregates for the collision are selected at random with a uniform distribution (Fig. B.1). In the linear kernel (Fig. B.2), the aggregates are selected weighted by the sum of the masses of the aggregates to collide, with a higher probability of collision for bigger aggregates. Analytical solutions (Tanaka & Nakazawa 1994; Silk & Takahashi 1979) for constant and linear kernels can be written as $$F\left(n\right)=N_{\text{agg,0}{g}^2}{\left(1-g\right)}^{n-1},$$(B.1) $$F\left(n\right)=N_{\text{agg,0}}\frac{n^{n-1}}{n!}g{\left(1-g\right)}^{n-1}{e}^{-n\left(1-g\right)}.$$(B.2)

Following the notation by Tanaka & Nakazawa (1994) g is the ratio of the total number of aggregates at a given time t (or n_coll) to its initial value (N_agg,0). g can also be expressed as a function of a normalized time τ = N_agg,0t, g = 1/(1 – τ/2) and e^–τ, for the constant and linear kernels, respectively. In this work, we take the time as t = n_collt₀, where t₀ is the average time between collisions, determined by the mean free path and the mean collision velocity, and assume constant t₀. In reality, t₀ will change dynamically with the dust and velocity distributions, depending on the specific astrophysical environment, resulting in various sequences of collision times. Here we present results for n, the number of particles in the aggregates rather than their mass m, remembering that this means a normalized mass because each monomer has the same mass.

Our results for F(n) are shown for the constant kernel in Fig. B.1 (a) and for the linear kernel in Fig. B.2(a). In our code, the number of aggregates is finite, starting with N_agg,0 monomers. Therefore, the distribution for large n in Fig. B.2(a) departs from the analytical one since full agglomeration occurs when the number of collisions approaches N_agg,0, as expected since a single aggregate will form, including all N_agg,0 monomers. If needed, this can be simply avoided starting with N_agg,0 larger than the desired relevant number of collisions, n_coll. Our results for F(n) n² are shown for the constant kernel in Fig. B.1(b) and for the linear kernel in Fig. B.2(b). We show excellent agreement with the analytical solutions for both kernels, including an accurate description of the distribution well beyond its maximum, with a sharp decrease for large aggregates F(n) n². These results can be compared with similar ones in Ormel et al. (2007), where F(m) m² was shown. In our case, a departure from the analytical results for n approaching N_agg,0 is observed, as explained above for panel (a).

Fig. B.1 Coagulation test with a constant kernel for (a) F(n) and (b) F(n)n2. Color lines show the analytic solution (Eq. (B.1)), and the dashed line with slope 1 is to guide the eye.

Fig. B.2 Coagulation test with a linear kernel for (a) F(n) and (b) F(n)n2. Color lines show the analytic solution (Eq. (B.2)), and the dashed line with slope 1 is to guide the eye.

Appendix C Mass distribution in different porosity ranges

In Sect. 4.4.1, different porosity ranges were analyzed. Fig. 6 presents the mass distribution, in logarithmic bins, for the three intervals of porosity P, PC, and C for three different values of ncoii. Alternatively, Fig. C.1 shows the mass distribution for each mentioned porosity interval and their differences when increasing n_coll. The differences in each porosity regime as n_coll grows show up clearly: For P case (Fig. a) none of the aggregates grows beyond the largest mass in the initial distribution, while in cases PC and C, Figs. (b) and (c), there is growth beyond that limit, with growth being larger for the larger number of collisions, as expected. In C case, Fig. C.1(c), there are aggregates with ϕ > 0.5 resulting from collisional processes, together with all monomers and dimers.

Fig. C.1 Aggregate mass distribution, F (normalized to unit area), for: (a) P, (b) PC, and (c) C porosity regimens. n indicates the number of monomers in a given aggregate, which is proportional to aggregate mass. Colors denote the distribution after ncoll collisions, with ncoll=0 indicating the initial distribution. Gray lines show fits of the form nα, with the dotted line for the low-mass region and the solid line for the high-mass region.

References

All Tables

All Figures

	Fig. 1 Snapshots at 50 µs showing the collision outcome – SP (sticking by penetration), TF (two fragments), or TD (total destruction) – between aggregates with µ = 10 and: (a) ϕ = 0.15, υ = 10 m s−1, (b) ϕ = 0.15, υ = 20 m s−1, and (c) ϕ = 0.40, υ = 50 m s−1.
In the text

	Fig. 2 Flow chart diagram of the Monte Carlo simulation.
In the text

	Fig. 3 Mass distribution (normalized by area). Colors denote the initial distribution and final distribution after 5, 10, and 100 collisions. Gray lines show fits for the exponent α (see text).
In the text

	Fig. 4 Size distribution (normalized by area). Colors denote the initial distribution and the final distribution after 5, 10, and 100 collisions. Gray lines show fits for the exponent β (see text).
In the text

	Fig. 5 Porosity distribution (normalized by area). Colors denote the initial distribution and the final distribution after 5, 10, and 100 collisions. (a) Complete distributions. (b) Only aggregates with n > 2.
In the text

	Fig. 6 Mass distribution (normalized by area) for ncoll: (a) 5, (b) 10, and (c) 100. Colors denote ranges of the porosity of the aggregates. Gray lines show fits for the exponent α (see text).
In the text

	Fig. 7 Size distribution (normalized by area) for ncoll: (a) 5, (b) 10, and (c) 100. Colors denote ranges of the porosity of the aggregates. Gray lines show fits for the exponent β (see text).
In the text

	Fig. 8 Size distribution (normalized by area) for ncoll = 100 at different Cconv. (a) Complete range of R. (b) 100 µm < R < 1 50 µm.
In the text

	Fig. A.1 Relation between initial conditions and final outcomes.
In the text

	Fig. B.1 Coagulation test with a constant kernel for (a) F(n) and (b) F(n)n2. Color lines show the analytic solution (Eq. (B.1)), and the dashed line with slope 1 is to guide the eye.
In the text

	Fig. B.2 Coagulation test with a linear kernel for (a) F(n) and (b) F(n)n2. Color lines show the analytic solution (Eq. (B.2)), and the dashed line with slope 1 is to guide the eye.
In the text

Fig. C.1 Aggregate mass distribution, F (normalized to unit area), for: (a) P, (b) PC, and (c) C porosity regimens. n indicates the number of monomers in a given aggregate, which is proportional to aggregate mass. Colors denote the distribution after ncoll collisions, with ncoll=0 indicating the initial distribution. Gray lines show fits of the form nα, with the dotted line for the low-mass region and the solid line for the high-mass region.

In the text