Issue 
A&A
Volume 536, December 2011



Article Number  A104  
Number of page(s)  11  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201117645  
Published online  19 December 2011 
Collisions of inhomogeneous preplanetesimals
^{1} Institut für Astronomie und Astrophysik, Abteilung Computational Physics, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: ralf.j.geretshauser@unituebingen.de
^{2} Physikalisches Institut, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany
Received: 7 July 2011
Accepted: 13 October 2011
Context. In the framework of the coagulation scenario, kilometresized planetesimals form by subsequent collisions of preplanetesimals of sizes from centimetre to hundreds of metres. Preplanetesimals are fluffy, porous dust aggregates, which are inhomogeneous owing to their collisional history. Planetesimal growth can be prevented by catastrophic disruption in preplanetesimal collisions above the destruction velocity threshold.
Aims. We assess whether the inhomogeneity created by subsequent collisions has a significant influence on the stability of preplanetesimal material to withstand catastrophic disruption. We wish to develop a model that is explicitly able to resolve any inhomogeneous structures. The input parameters of this model must be easily accessible from laboratory measurements.
Methods. We develop an inhomogeneity model based on the density distribution of dust aggregates, which is assumed to be a Gaussian distribution with a welldefined standard deviation. As a second input parameter, we consider the typical size of an inhomogeneous clump. For the simulation of the dust aggregates, we utilise a smoothed particle hydrodynamics (SPH) code with extensions for modelling porous solid bodies. The porosity model was previously calibrated for the simulation of SiO_{2} dust, which commonly serves as an analogue for preplanetesimal material. The inhomogeneity is imposed as an initial condition on the SPH particle distribution. We carry out collisions of centimetresized dust aggregates of intermediate porosity. We vary the standard deviation of the inhomogeneous distribution at fixed typical clump size. The collision outcome is categorised according to the fourpopulation model.
Results. We show that inhomogeneous preplanetesimals are more prone to destruction than homogeneous aggregates. Even slight inhomogeneities can lower the threshold for catastrophic disruption. For a fixed collision velocity, the sizes of the fragments decrease with increasing inhomogeneity.
Conclusions. Preplanetesimals with an active collisional history tend to be weaker. This is a possible obstacle to collisional growth and needs to be taken into account in future studies of the coagulation scenario.
Key words: hydrodynamics / methods: numerical / planets and satellites: dynamical evolution and stability / protoplanetary disks / planets and satellites: formation
© ESO, 2011
1. Introduction
Terrestrial planets are thought to form by core accretion in protoplanetary accretion discs (e.g. Lissauer & Stewart 1993). The discs consist of both hydrogen and helium gas and a small fraction of dust, which in the early stages exists as micron sized dust grains (e.g. Dullemond et al. 2007). Interaction between the gas and dust induces relative motions between the dust grains and causes collisions among them (e.g. Weidenschilling & Cuzzi 1993). Mutual sticking by van der Waals forces occurs and by this process millimetre to centimetresized, highly porous preplanetesimals form (e.g. Blum & Wurm 2008). Additional growth by means of collision is endangered by several obstacles: rebound, material loss by accretion or photoevaporation, and fragmentation. Zsom et al. (2010) found that the growth of preplanetesimals to kilometresized planetesimals may be halted by mutual compaction and bouncing, which leads finally to a population of dust and pebblesized objects in the protoplanetary disc. Weidenschilling (1977) found that metresized bodies have the highest drift speed towards the host star, where the objects get accreted or photoevaporated. For these objects, the drift timescale is much shorter than the time needed for a planet to form. One of the most serious obstacles to planetesimal formation is the fragmentation barrier: with increasing aggregate size the mutual collision velocities also increase and catastrophic disruption might become more common than mass gain (e.g. Brauer et al. 2008).
If by some mechanism a sufficient population of kilometresized objects is formed, their selfgravity acts as an accretion mechanism and ensures further growth to planets. In this gravity dominated regime, kilometresized objects and yet larger ones are referred to as planetesimals.
To investigate planetesimal formation, global dust coagulation models (e.g. Dullemond & Dominik 2005; Zsom et al. 2010; Birnstiel et al. 2010) numerically compute the evolution of dust aggregates in a protoplanetary disc starting from micronsized dust grains. These models require detailed collision statistics evaluated for the outcome of twobody preplanetesimal collisions as functions of important parameters such as mass ratio, porosity, impact parameter, and relative velocity of the collision partners. In the past decade, many experimental studies have been carried out to collect these data (e.g. Blum & Wurm 2008; Güttler et al. 2010, for summaries). However, laboratory setups become difficult and infeasible for centimetresized preplanetesimals and larger ones since they have to be carried out under protoplanetary disc conditions, i.e. in a vacuum and under microgravity. In this regime, collision data has to be acquired by numerical simulations.
In their pioneering work, Benz & Asphaug (1994, 1995) developed a solid body smoothed particle hydrodynamics (SPH) implementation for brittle material based on the dynamical fracture model by Grady & Kipp (1980), which evolves the propagation of cracks starting from an initial Weibull (1939) distribution of flaws. On these grounds, Benz (2000) investigated collisions between rocky nonporous preplanetesimals in their typical collision velocity regime from 5 m s^{1} to 40 m s^{1} and found only disruption. However, experiments (e.g. Blum & Wurm 2008) indicate that preplanetesimals are fluffy, porous objects rather than solid rocks. For this reason, Sirono (2004) developed an SPH implementation of a porosity model based on porositydependent material quantities such as bulk and shear moduli as well as compressive, shear, and tensile strengths. In his simulations of porous ice aggregates, he found that owing to their porous structure the kinetic energy can be dissipated effectively, permitting collisional growth. Sirono’s approach included a damage model based again on the dynamical fracture model of Grady & Kipp (1980). In addition, a damage restoration model tries to capture damage healing by compression. However, Schäfer et al. (2007) found that Sirono’s damage model is not applicable to their SPH simulations of porous ice because it includes damage by compression. In contrast, materials such as porous ice and porous dust form new molecular bondings when they are compressed. Thus, their strengths increase but do not decrease as predicted by the Sirono damage model. On the basis of the implementation of Benz & Asphaug (1994, 1995), Jutzi et al. (2008) combined the damage model of Grady & Kipp (1980) and the p − α model of Herrmann (1969) to simulate porous brittle media. The model was calibrated with laboratory high velocity impact experiments (Jutzi et al. 2009). However, the model by Jutzi et al. (2008) is also not applicable to porous protoplanetary dust because it includes damage by compression and has no suitable model of damage healing.
To simulate SiO_{2} dust, which is used as protoplanetary dust analogue in many experimental studies, we developed a suitable porosity model based on the ideas of Sirono (2004), that was implemented into a solid body SPH code to numerically investigate preplanetesimal collisions (Geretshauser et al. 2010). With the aid of laboratory benchmark experiments (see also Güttler et al. 2009), the porosity model was calibrated and tested for the simulation of SiO_{2} dust. In particular, the simulation reproduced the compaction, bouncing, and fragmentation behaviour of homogeneous dust aggregates of various porosities in a quantitatively correct way. However, this approach does not include an explicit damage model.
In this paper, we present an approach to explicitly resolve inhomogeneous structures of dust aggregates. Since an inhomogeneous density distribution is one of the main sources of fracture in porous aggregates, in this work we approximate the damage evolution by the evolution of density. The inhomogeneity is described by a Gaussian density distribution and the typical size of a clump. The aggregate inhomogeneity of SiO_{2} dust can easily be measured in the laboratory by Xray tomography (see Güttler et al. 2009, Fig. 10). In contrast, the Grady & Kipp (1980) damage model relies on the correct parameters for the initial Weibull distribution of flaws, which are difficult to determine in laboratory experiments. Our approach of resolving inhomogeneous structures is closely connected to the porosity model, where the strength quantities are porosity dependent. Therefore, any fluctuation in porosity induces a fluctuation in strength over the aggregate, which is utilised to model the evolution of the material under stress.
To classify the collisional outcome of our simulations involving inhomogeneous aggregates, we utilise the fourpopulation model. This classification scheme was introduced by Geretshauseret al. (2011) to enable a precise mapping of all collisional outcomes involving any combination of sticking, bouncing, and fragmentation processes. Four fragment populations are distinguished according to their masses: both the largest and the second largest fragment as well as a set of fragments whose mass distribution can be approximated by a power law. To take into account the resolution limit, the subresolution population contains all fragments that consist of a single SPH particle.
The outline of this paper is the following. In Sect. 2, we describe the SPH method with extensions for the simulation of solid bodies and the applied porosity model. Section 3 firstly describes a previous study of a damage model, then we present our approach based on the inhomogeneity of dust aggregates and its implementation. The results of collisions with inhomogeneous aggregates are discussed in Sect. 4. We summarise our findings in Sect. 5, where we also discuss the results and provide an outlook on future work.
2. Solid body SPH and porosity model
For our simulations, we apply the smoothed particle hydrodynamics (SPH) numerical method originally developed by Lucy (1977) and Gingold & Monaghan (1977) for astrophysical flows. In subsequent studies, SPH was extended to model the elastic and plastic deformations of solid materials (e.g. Libersky & Petschek 1991; Benz & Asphaug 1994; Randles & Libersky 1996). Our parallel SPH code parasph (Hipp & Rosenstiel 2004; Schäfer 2005; Schäfer et al. 2007; Geretshauser et al. 2010) is based on these concepts. In the SPH scheme, the time evolution of a solid body is simulated by means of the Lagrangian equations of continuum dynamics which represent the continuum and momentum equations, respectively. The quantities ρ, v_{α}, and σ_{αβ} have their usual meanings, i.e. density, velocity, and the stress tensor. The Greek indices denote the spatial components and the Einstein sum convention applies to them. Since our equation of state (EOS) is energyindependent, we do not solve the energy equation, hence omit to state it here. Within the framework of SPH, the continuous quantities are discretised into mass packages, which we refer to as “particles”. These represent the sampling points of the method and interact with each other within a limited spatial range, called the smoothing length h. The influence of a particle b on a particle a depends on the particle distance x^{a} − x^{b} and is weighted by the kernel function W(x^{a} − x^{b},h). This is usually a spherically symmetric function with compact support, which is differentiable at least to first order. We utilise the cubic spline kernel of Monaghan & Lattanzio (1985).
The particles of the SPH scheme move according to the Lagrangian form of the equations of continuum mechanics. Thus, they represent a natural frame of reference for any deformation of the simulated solid body. In the context of fragmentation, two fragments are considered to be separate once their constituent subsets of particles no longer interact, i.e., when the fragments are separated by more than 2 h. The fragments may be spread out over an unlimited computational domain. This is why computations are carried out at the particle positions.
The stress tensor in the momentum equation can be divided into the term accounting for pure hydrostatic pressure p and the term representing pure shear, which is given by the traceless deviatoric stress tensor S_{αβ}. Hence, (3)In contrast to viscous fluids, p and S_{αβ} do not depend on the velocity gradient but on the deformation of a solid body given by the strain tensor (4)where x′ are the coordinates of the deformed body. In particular, the relations that hold are where the quantity d denotes the dimension and K and μ are the bulk modulus and shear modulus, respectively, for elastic deformation. To obey the principle of frame invariance, rotation terms have to be added to the deviatoric stress tensor. For this, we apply the Jaumann rate form, which was previously used by Benz & Asphaug (1994) and Schäfer et al. (2007), and described by Geretshauser et al. (2010) in more detail.
The above relations (Eqs. (5) and (6)) describe the time evolution of a solid body in the elastic regime, i.e., for reversible deformations. With respect to preplanetesimal collisions, we are concerned with highly porous material that undergoes elastic and plastic deformations. Both aspects are described by means of the porosity model of Geretshauser et al. (2010), which is based on the ideas of Sirono (2004) and can be categorised as a plasticitybased porosity model.
In this approach, the quantities for the elastic regime (bulk and shear moduli), as well as the strength quantities (compressive, shear, and tensile strength), depend on the filling factor, which accounts for the subresolution porosity. The filling factor φ is related to the density ρ and the porosity Φ by (7)where ρ_{s} is the density of the matrix material. We use a constant value of ρ_{s} = 2000 kg m^{3} for our SiO_{2} material (e.g. Blum & Schräpler 2004). In particular, the elastic quantities are related to the filling factor by a power law (8)where γ = 4, K_{0} is the bulk modulus of an uncompressed random ballistic deposition (RBD) dust sample with φ_{RBD} = 0.15 (e.g. Blum & Schräpler 2004), such that K_{0} = K(φ_{RBD}) = 4.5 kPa.
For the compressive, tensile, and shear strength we adopt some empirical relations (see Güttler et al. 2009; Geretshauser et al. 2010; and Geretshauser et al. 2011, for details). These relations allow us to reproduce the compaction, bouncing, and fragmentation behaviour to high accuracy in a velocity range of the order of 0.1 to 10 m s^{1} (Geretshauser et al. 2010). Geretshauser et al. (2011) showed that all collision types found in laboratory experiments can be reproduced. In the same reference, the laboratory velocity thresholds for the transitions between sticking, bouncing, and fragmentation behaviour resembled the simulation results for collisions of medium porosity aggregates with up to 27.5 m s^{1}. Therefore, we expect these relations to be valid for collisions below the sound speed of the dust material, which is ~ 30 m s^{1} (Güttler et al. 2009). For the compressive strength, (9)where φ_{min} + ε < φ < φ_{max} and ε = 0.005. The quantities φ_{min} = 0.12 and φ_{max} = 0.58 are the minimum and maximum filling factors, respectively, in the compressive strength relation. Both values can be exceeded by the material. The power of the compressive strength relation is ln(10) times the parameter Δ with Δ = 0.58. The quantity p_{m} = 260 Pa accounts for the mean pressure of the compressive strength relation. According to Geretshauser et al. (2010, 2011), these parameters are valid for impacts below ~30 m s^{1}. The influence of Δ and p_{m} was studied by Geretshauser et al. (2010). For φ ≤ φ_{min} + ε, the compressive strength relation is continuously extended by the constant function Σ(φ) = Σ(φ_{min} + ε) and for φ_{max} ≤ φ, we set Σ(φ) = ∞. The tensile strength is given by (10)and the shear strength is the geometric mean of the compressive and tensile strength (11)We note that the three strength quantities are fillingfactor dependent. This is because as the filling factor increases the monomers of the dust material are packed more closely. Thus, each monomer establishes more bonds, which increases the strength. The fillingfactor dependence of the strength quantities is exploited by the approach presented in this paper. Together the bulk modulus, the compressive strength, and the shear strength yield the full equation of state for pure hydrostatic compression or tension (12)where , and and are critical filling factors. The value of marks the transition between elastic and plastic compression and defines the transition between elastic and plastic tension. The filling factors inbetween the critical values represent the elastic regime. The quantity is the reference filling factor that represents the filling factor of the porous material at vanishing external pressure. Once the material undergoes plastic deformation, the reference filling factor is reset using the middle relation of Eq. (12). Details of this pressure reduction process can be found in Geretshauser et al. (2010). In the elastic regime (middle line of Eq. (12)), the pressure varies with density according to the elastic path given by Eq. (8). If is exceeded, the pressure follows the much shallower compressive strength path Σ(φ) in the positive φdirection, whereas if is exceeded, the pressure follows the shallower tensile strength path T(φ) in the negative φdirection (cf. Fig. 1 in Geretshauser et al. 2010). In both cases, the absolute value of the pressure, which is computed according to the elastic equation of state, is reduced to the compressive or tensile strength, respectively. This accounts for the stress relieved by plastic deformation. Details of the implementation are described by Geretshauser et al. (2010).
For the plastic deformation by pure shear, the deviatoric stress tensor S_{αβ} is reduced according to the von Mises plasticity given by (e.g. Benz & Asphaug 1994) (13)The function f is defined by . An alternative approach (e.g. Wilkins 1964; Libersky et al. 1993) uses the expression but we expect this to lead to differences from our adopted approach that are marginal. In this relation, J_{2} = 0.5 S_{αβ}S_{αβ} and Y = Y(φ) are the second irreducible invariant of the deviatoric stress tensor and the shear strength, respectively.
3. Resolving inhomogeneous structures
3.1. Previous approach
In contrast to ductile media, brittle materials, such as basalt, granite, or porous pumice, do not rupture by plastic flow. This is because the material is not completely homogeneous but contains little flaws. These are little defects in the medium. With increasing strain, cracks develop originating from these flaws and start to pervade the solid body. In brittle media, stress is relieved by developing cracks.
To model damage, Grady & Kipp (1980) introduce a scalar parameter 0 ≤ D ≤ 1, where D = 0 represents undamaged and D = 1 fully disintegrated material. The stress tensor (Eq. (3)) is modified according to (14)Thus, damaged material may feel less stress than undamaged material. The local damage D is defined as the ratio of the volume defined by the growing crack to the volume in which the crack is growing (see Grady & Kipp 1980). The local time evolution of the damage D is given by (15)where n_{act} accounts for crack accumulation and denotes the number of activated flaws, and R_{s} is the radius of a circumscribing sphere in which the crack evolves. It accounts for the maximum size of a damaged area. The quantity c_{g} is the crack growth velocity. The concept of flaw activation originates from the idea that cracks do not start to grow for any strain applied but are activated for some strain threshold ϵ_{act}. The number of flaws n per unit volume with ϵ_{act} is given mostly by a power law (16)This distribution of flaws in a brittle material was proposed by Weibull (1939). It is based on two material parameters k_{wb} and m_{wb}, where k_{wb} is the number of flaws per unit volume. The flaw distribution according to Eq. (16) is set as an initial condition for the material before the simulation.
A disadvantage is that the material parameters k_{wb} and m_{wb} are rarely available because they are difficult to measure. Unfortunately, small variations in k_{wb} and m_{wb} lead to large differences in the activation thresholds and simulation results (e.g. Schäfer 2005; Jutzi et al. 2009).
Because the damage model by Grady & Kipp (1980) is only suitable for brittle material, it is not applicable to preplanetesimal material represented by SiO_{2} dust, which has properties that are between those of a ductile and a brittle material. Therefore, it is desirable to develop a simple damage approximation consistent with the applied porosity model. To summarise, three ingredients are essential for a damage model: (1) a distribution of defects, (2) the stress reduction due to damage, and (3) the typical size of a damaged area.
3.2. The approach of resolving inhomogeneous structures
In contrast to the damage model presented in Sect. 3.1, our approach does not explicitly consider the evolution of defects, hence is not strictly a damage model. However, it is designed to behave as an approximation to a damage model in some aspects, which are described in this section. As a significant advantage, our approach relies on quantities that can be more easily measured by laboratory experiments than the parameters of the model by Grady & Kipp (1980).
The starting point for our approach is the inhomogeneous nature of the material. Macroscopic dust aggregates created by the RBD method are not completely homogeneous (Güttler et al. 2009). The filling factor instead was found to follow a Gaussian distribution around a median of φ_{μ} ~ 0.15 with φ_{σ} = 0.013.
According to the porosity model of Sect. 2, regions of lower filling factor also represent regions of weaker compressive, tensile, and shear strengths, whereas regions of higher filling factor are stronger. Therefore, as an analogue to the Weibull distribution we propose an initial distribution of the filling factor given by the Gaussian function (17)where n(φ) is the number density for a filling factor φ, φ_{σ} is the width of the distribution function, and φ_{μ} is the median filling factor. The pressure limits of compressive, tensile, and shear strengths Σ(φ), T(φ), and Y(φ), respectively, mark the transition from the elastic to the plastic hydrostatic regime. Therefore, they can be regarded as analogues to the activation threshold ϵ_{act}. Hence, we can associate with each φ analogues to the activation threshold, which accounts for the first criterion of a damage model (see end of Sect. 3.1) that is based on the distribution of defects.
Damaged regions in this simple inhomogeneity scheme are therefore represented by areas of low filling factor. Under constant tension, φ becomes increasingly smaller with time in these regions and, owing to the φdependence of either the tensile strength T(φ) or shear strength Y(φ), the strength decreases in these regions as in Eq. (14) with increasing D. This fulfils the second criterion: the stress reduction.
Instead of artificially introducing a damage parameter D and an activation threshold ϵ_{act}, the filling factor φ naturally takes over this twofold role. On the one hand, it represents a quantity which determines the activation thresholds for the plastic regime, and on the other, it also represents a damage parameter that reduces the strength quantities. Therefore, in the inhomogeneity approach, the time evolution of the approximated damage is given by the time evolution of φ or equivalently the density. The propagation speed of the defects is, intrinsically, the sound speed.
We introduce a typical size R_{c} of an inhomogeneous clump to comply with the third criterion. In terms of this parameter, an absolute length scale is introduced into the inhomogeneity approach.
A great advantage of our approach is that the material parameters φ_{σ}, φ_{μ}, and R_{c} can be determined in laboratory measurements, e.g. by Xray tomography (Güttler et al. 2009). However, no systematic study has been performed yet. A single measurement for a random ballistic deposition dust sample, which has the lowest filling factor (φ = 0.15) and the highest degree of homogeneity, has found a standard deviation of φ_{σ} ~ 0.013 (Fig. 10 in Güttler et al. 2009). There exist no data on either the typical clump size or clump size distribution. Therefore, as a first approximation, in Sect. 4 we study the effect of φ_{σ} on the outcome of preplanetesimal collisions for a fixed clump size R_{c}. The influence of the latter will be studied in future work.
3.3. Implementation issues
The inhomogeneity is imposed on the simulated aggregate as an initial condition. Initially, we generate a homogeneous aggregate by assigning a constant initial filling factor φ_{i}. As a second step, φ_{i} is modified according to a Gaussian distribution with a standard deviation φ_{σ}, which we will simply call either the “inhomogeneity” or “degree of inhomogeneity” in this paper. For the distribution, we randomly select a particle a and set its density to a new following the Gaussian. The same is assigned to all particles within a radius of typical clump, which is fixed to R_{c} = 9.75 mm, the equivalent of one smoothing length h. Together with a spatial distance of 2.6 mm between the SPH particles, one clump is resolved by ~290 particles, which is well above the required resolution (Geretshauser et al. 2010). Since there are no empirical data for the typical clump size yet, as a first approximation we assume that the clump size is constant. With the chosen size, the simulated spherical decimetre bodies can be regarded as being assembled from roughly spherical centimetre aggregates in sticking processes with different velocities similar to those described by Güttler et al. (2010) and Geretshauser et al. (2011). However, the effect of the typical clump size and resolution on the resulting fragment distribution will be studied in future work.
It is very likely that dust aggregates are compacted during their evolution in the protoplanetary disc (Zsom et al. 2010), in particular if the aggregates are assembled by smaller units in sticking processes. For this reason, we do not consider very fluffy aggregates with φ_{i} = 0.15 as produced by random ballistic deposition. We instead simulate spheres whose filling factor lies between the minimum (φ_{min} = 0.12) and maximum (φ_{max} = 0.58) filling factor of the compressive strength relation (Eq. (9)), i.e. φ_{i} = 0.35.
The mass of all individual SPH particles, which is a numerical parameter in the SPH equations, takes a constant value in the simulation. This value was computed to be consistent with the initial homogeneous aggregate (here φ_{i} = 0.35). Therefore, the modification of the density introduces a slight inconsistency in the SPH particle distribution. This means that the density value of an arbitrary particle differs slightly from the value computed from the particle number density using the standard SPH sum approach. Potentially, this could introduce pressure fluctuations causing spurious particle motion. However, in our algorithm the density is determined by solving the continuity equation and not by evaluating the number density of the particles (e.g. Monaghan 2005). As a consequence this inconsistency is marginal. Moreover, the time evolution of a inhomogeneous dust aggregate at rest was simulated and no spurious particle motion has been detected. In addition, the aggregate displayed no signs of instabilities within simulation times much longer than the collision timescale. Consequently, the presented approach is assumed to yield stable aggregates. Evaluating the generated SPH particle distribution yields the resulting filling factor distribution: the number density n(φ) of the volume fraction with filling factor φ is shown in Fig. 1 as a binned and cumulative plot for various φ_{σ}.
Fig. 1 Initial filling factor distributions of the target and projectile as binned (top) and cumulative diagram (bottom). The number density n(φ) of volumes with filling factor φ is plotted over φ. A larger degree of inhomogeneity is characterised by the increasing standard deviation φ_{σ} of the Gaussian (Eq. (17)) around the initial filling factor φ_{i} of the homogeneous aggregate. Here, we choose φ_{i} = 0.35. The values of φ_{σi} can be found in Table 1. 
Fig. 2 Interior and exterior view of the targets for different inhomogeneities with the filling factor colour coded. The inhomogeneity increases from the homogeneous target a) to the most inhomogeneous target f). In particular, the standard deviations are φ_{σ} = 0 a), φ_{σ} = 0.01 b), φ_{σ} = 0.02 c), φ_{σ} = 0.03 d), φ_{σ} = 0.04 e), and φ_{σ} = 0.05 f). All targets have a similar pattern that originates from considering interacting particle packages in the implementation. With increasing φ_{σ} the maximum and minimum filling factor of the different spots increase and decrease, respectively. The projectiles have a similar appearance. The result of collisions among these aggregates with 10 m s^{1} are depicted in Fig. 3. 
The inhomogeneous clumps with different filling factors are clearly visible in Fig. 2, where the density structure of the target is colour coded for different φ_{σ}. The packages are picked randomly but the random pattern is the same for all aggregates to ensure that the results are comparable. Owing to the increasing φ_{σ}, the regions with φ ≠ φ_{i} gain higher maximum and lower minimum filling factors. Consequently, the inhomogeneity pattern for φ_{σ1} and φ_{σ2} (see Table 1) is barely visible in Fig. 2 (first row, second and third). In contrast, for φ_{σ5} the difference between the extreme values of the filling factor create a clearly visible pattern with high contrast.
The standard deviation of the Gaussian filling factor curve φ_{σ} is supplied to the distribution routine as an input parameter. The range of φ_{σ} begins slightly below the empirical value for a homogeneous random ballistic deposition aggregate (φ_{σ} = 0.013). There is no empirical data on the inhomogeneity of aggregates resulting from subsequent collisions. However, owing to the collisional history of an aggregate with macroscopic preplanetesimals it is very likely that realistic aggregates have a higher degree of inhomogeneity. Therefore, we increase φ_{σ} in steps of 0.01 up to φ_{σ} = 0.050. The real standard deviation computed from the full width at half maximum (FWHM) is smaller than the φ_{σ} supplied to the SPH particle initialisation algorithm (algorithmic value). The FWHM values are listed in Table 1. For the sake of simplicity, we use the algorithmic value this paper.
Selection of the standard deviations φ_{σ} as they are used in the inhomogeneity algorithm.
4. Simulation results
We now perform test simulations using our inhomogeneity approach. The initial setup is the same as in most simulations of Geretshauser et al. (2011). It consists of two colliding spheres with φ_{i} = 0.35. The radii of the target and projectile are r_{t} = 10 cm and r_{p} = 6 cm, and they consist of 238 238 and 51 477 SPH particles, respectively, with masses of 1.23 × 10^{5} kg each. The SPH particles are set on a cubic grid separated by a distance of 2.6 mm. We consider two impact velocities of v_{0} = 10 m s^{1} and v_{0} = 12.5 m s^{1}. For homogeneous aggregates, the smaller velocity leads to the sticking of the projectile to the target (see Fig. 3a). In contrast, the higher velocity results in fragmentation (see also Geretshauser et al. 2011). For these two collision velocities, the inhomogeneity of both the target and the projectile is defined by the standard deviation of the Gaussian φ_{σ}, which is varied from 0 to 0.05 yielding the distributions depicted in Fig. 1. The resulting targets are shown in Fig. 2. The projectiles display a similar density pattern.
Fig. 3 Outcome of a collision between a target and projectile with the same φ_{σ} for different inhomogeneities. In all cases, the target and projectile radii are r_{t} = 10 cm and r_{p} = 6 cm, respectively, and the collision velocity is v_{0} = 10 m s^{1}. Analogous to the initial targets in Fig. 2, the standard deviations for the Gaussian of the initial targets are φ_{σ} = 0 a), φ_{σ} = 0.01 b), φ_{σ} = 0.02 c), φ_{σ} = 0.03 d), φ_{σ} = 0.04 e), and φ_{σ} = 0.05 f). The collision outcome is shown in the impact direction. In the homogeneous case a) and for small inhomogeneities b), the target stays intact and forms one massive object with the projectile. For φ_{σ} ≥ 0.02, the target fragments c)–f). The fragment sizes decrease with increasing φ_{σ} and at the same time the number of fragments increases. 
For v_{0} = 10 m s^{1}, the final fragment distribution is shown in Fig. 3 from the impact direction. For homogeneous aggregates (a) and the small standard deviation (φ_{σ} = 0.01, b), the target remains intact and the target and projectile form one large aggregate. For φ_{σ} = 0.01, some small fragments are visible. For φ_{σ} ≥ 0.02, the target fragments and with greater inhomogeneity the largest fragments of the distribution decrease in size.
The first impression from our visual inspection of the collision products is confirmed by the more detailed analysis. For this purpose, we utilise the fourpopulation model (Geretshauser et al. 2011). According to this classification scheme of collision outcome, four classes of fragments are distinguished according to their mass: the largest fragment (index “1”) and the second largest fragment (index “2”) are characterised by single values of mass m, kinetic energy E, and other parameters. To account for the resolution limit, the subresolution population (index “sr”) is introduced, which contains all fragments consisting of a single SPH particle (mass m_{SPH}). This class is mainly described by global quantities. Fragments with masses m_{SPH} < m_{frag} < m_{2} form a class whose cumulative mass distribution is described by a powerlaw (index “pw”). This class is characterised by global quantities and distribution functions.
The evolution of the masses of each population with φ_{σ} is shown in Fig. 4 as a fraction of the total mass. The upper and the lower plots display the results for v_{0} = 10 m s^{1} and v_{0} = 12.5 m s^{1}, respectively. The exact values can be found in Table 2. For the lower collision velocity and φ_{σ} ≤ 0.01, nearly all mass is stored in the largest fragment m_{1}. For φ_{σ} > 0.01, a second largest fragment and a powerlaw population appears. With increasing inhomogeneity m_{1} rapidly decreases, while the mass of the second largest fragment m_{2} only slightly decreases. The masses m_{1} and m_{2} become comparable in the end. For high inhomogeneity values, most of the mass is stored in the powerlaw population m_{pw}, as we discuss below. However, for high inhomogeneities m_{pw} seems to remain fairly constant, while the mass of the subresolution population m_{sr} increases. For the higher collision value (bottom plot of Fig. 4), the evolution is similar but starts with m_{pw} ≠ 0. The mass m_{1} constantly decreases, but not as rapidly as in the low velocity case. Moreover, m_{2} at first increases, then slightly decreases until the largest and second largest fragment become comparable. The mass of the powerlaw population rapidly increases for 0.01 ≤ φ_{σ} ≤ 0.04. For larger φ_{σ}, the powerlaw population only slightly increases, whereas m_{sr} increases at a higher rate. The mass evolution with increasing inhomogeneity can be summarised as follows: (1) inhomogeneity leads to fragmentation at collision velocities, at which homogeneous aggregates do not fragment, i.e. inhomogeneous aggregates are weaker; (2) for the two investigated velocities a larger φ_{σ} leads to a decrease in the mass of the largest and second largest fragments; (3) a higher degree of inhomogeneity causes an increase in the masses of both the powerlaw and subresolution population.
We turn to the analysis of the distribution of the residual total energy after the collision, which is shown in Fig. 5 as a fraction of the initial kinetic energy. The residual total energy is the sum of the translational energy and the energy stored in the internal degrees of freedom (e.g. rotation, vibration) of the respective populations. The tabulated values can be found in Table 2. In Fig. 5, the low velocity case with v_{0} = 10 m s^{1} is shown in the upper plot and the high velocity case with v_{0} = 12.5 m s^{1} in the lower plot. For the former, a large amount of energy E_{1} is stored in the largest fragment, but a considerable amount E_{sr} is also stored in the subresolution population. As E_{1} decreases, the contribution of the second largest fragment at first increases, then decreases slightly with larger φ_{σ}. As for the mass, the energy contribution of the powerlaw population E_{pw} increases with a higher degree of inhomogeneity. Remarkably, the total energy fraction of these three populations remains almost constant. In contrast, the energy contribution of the subresolution population E_{sr} increases more strongly than its mass contribution in Fig. 4. The overall residual energy increases with a higher degree of inhomogeneity and this energy excess is stored mostly in the subresolution population. For the high velocity case (bottom plot of Fig. 5), this is even more evident: E_{1}, E_{2}, and E_{pw} behave similar to the low velocity case but with E_{2} ≠ 0 and E_{pw} ≠ 0 for φ_{σ} = 0. For a higher degree of inhomogeneity, less energy can again be dissipated and the total residual energy of the largest three populations remains nearly constant for increasing φ_{σ}. As a consequence, the residual energy excess is mainly stored in the subresolution population for large values of φ_{σ}. The conclusions of our energy analysis are: (1) with increasing impact velocity less energy can be dissipated by the system and the residual energy increases as previously found by Geretshauser et al. (2011); (2) with increasing inhomogeneity, the effectiveness of the energy dissipation decreases for the two investigated velocities; (3) the energy contribution of all populations show a similar dependence on φ_{σ} to their mass contributions: the fractions of energy stored in the largest and second largest fragment follow a decreasing trend with increasing inhomogeneity. The contribution of powerlaw and subresolution populations increase; (4) while the overall energy contribution of the largest fragment populations remains nearly constant with φ_{σ}, the fraction stored in the subresolution population drastically increases, such that the energy that cannot be dissipated owing to the increasing inhomogeneity is stored in rapidly moving single SPH particles.
Fig. 4 The total masses of the four populations with inhomogeneities measured by φ_{σ}. The collision velocities are v_{0} = 10 m s^{1} (top) and v_{0} = 12.5 m s^{1} (bottom). The masses of the largest fragment m_{1}, of the second largest fragment m_{2}, of the powerlaw population m_{pw}, and of the subresolution population m_{sr} are stacked and normalised by the total mass of the system m_{tot}. In both cases, m_{1} and m_{2} decrease with increasing φ_{σ}. Conversely, m_{pw} and m_{sr} increase. For the homogeneous case (φ_{σ} = 0), sticking is found for the low velocity collision, whereas for v_{0} = 12.5 m s^{1}m_{2} ≠ 0 and m_{pw} ≠ 0, which indicates fragmentation. 
Fig. 5 The energy of each population for increasing inhomogeneity parameter φ_{σ} and collision velocity v_{0} = 10 m s^{1} (top) and v_{0} = 12.5 m s^{1} (bottom). The energy of each population is the sum of translational, rotational and vibrational energy and normalised by the initial kinetic energy E_{init}. The two latter energy contributions are negligible for the headon collisions presented here. Owing to the strong correlation with mass, this diagram follows a similar trend to Fig. 4. The energy contributions are divided up into E_{1}, E_{2}, E_{pw}, and E_{sr} for the largest and second largest fragment, the powerlaw and subresolution population, respectively. In both velocity cases, both E_{1} and E_{2} decrease with increasing φ_{σ}. In contrast, E_{pw} and E_{sr} increase. In particular, the high fraction of kinetic energy stored in the subresolution population is remarkable. 
In comparison to the initial filling factor distribution of Fig. 1, the final filling factor distribution of the cumulated mass is shifted towards higher filling factors. Figure 6 reveals that in most cases more than 85% of the fragment mass is stored in filling factors ≳ 42% (except for v_{0} = 12.5 m s^{1} and φ_{σ} = 0.050). With the same exception, the distributions look quite similar in both velocity cases. There is a slight tendency for there to be more mass at lower filling factors in the collision with v_{0} = 12.5m s^{1}. From Fig. 6, it can also be seen that the broader the initial filling factor distribution, the broader the final filling factor distribution. For the v_{0} = 10.0 m s^{1} case, the two most homogeneous aggregates hardly produce any fragments. This causes a jump in the curves without any intermediate steps. The distribution for the highest collision velocity (v_{0} = 12.5 m s^{1}) and the highest degree of inhomogeneity (φ_{σ} = 0.050) is particularly interesting since ~30% of the fragment mass possess filling factors <0.35, which is a larger fraction than for the initial aggregate. This implies that rarefaction occurs during the collision. On the microscopic scale, this is equivalent to monomer reordering and stretching of monomer chains. We note that a tiny fraction of particles acquire overdensities (φ > 1) in the collision. In the most violent case, this applies to single SPH particles, which in total represent 0.03% of the total mass of the system, which is negligible. These overdensities indicate that the monomers are breaking up on the microscopic scale.
We have so far considered all members of the fourpopulation model. We now constrain our analysis to the powerlaw population. Its mass distribution is shown in the cumulative plot of Fig. 7 again for v_{0} = 10 m s^{1} (top) and v_{0} = 12.5 m s^{1} (bottom). The figure shows the ratio of the cumulative mass m_{cum} to the fragment mass m_{frag}, which are both normalised by the total mass of the powerlaw population m_{pw} given in Table 2. The cumulative distribution can be accurately (e.g. Güttler et al. 2010) described by a power law (18)where n(m) dm is the number density of fragments of mass m, which is normalised by the total mass of the powerlaw population m_{pw}. The quantity μ_{pw} is the normalised mass of the most massive member of the powerlaw population and κ is the powerlaw index or fragmentation parameter. We fitted the cumulative mass distributions of Fig. 7 for all fragment masses with Eq. (18). Finding an explanation for the kink in the cumulative mass curves is the subject of ongoing research. This deviation from a power law might be either a physical or a resolution effect (see also Geretshauser et al. 2011). The results are listed in Table 3. Figure 8 shows that the powerlaw index κ decreases slightly as the degree of inhomogeneity increases for both collision velocities, which leads to increasingly shallower slopes in Fig. 7 owing to the production of smaller fragments with increasing φ_{σ}. In addition, the κ values for the higher impact velocity (v_{0} = 12.5 m s^{1}) are systematically below the κ values of the lower velocity (v_{0} = 10.0 m s^{1}), except for small values of φ_{σ} for which there is only a small number of fragments (see Table 2) and hence the quality of the statistics is insufficient. This finding indicates that a larger fraction of less massive fragments are produced at higher impact velocities, which is expected.
Fig. 6 Cumulative mass of the fragment masses over their average filling factor φ. The mass is normalised by the total fragment mass m_{tot}. The collision velocities are v_{0} = 10 m s^{1} (top) and v_{0} = 12.4 m s^{1} (bottom). Compared to the initial distribution of filling factors (Fig. 1), the curves are shifted significantly towards higher filling factors. This reflects the compression taking place during the collision. 
Fig. 7 Cumulative mass distribution of a powerlaw population for different inhomogeneity parameters φ_{σ} and collision velocities v_{0} = 10 m s^{1} (top) v_{0} = 12.5 m s^{1} (bottom). In both velocity cases both the mass of the largest member of the powerlaw population and the powerlaw slope decrease. This indicates the fragmentation to smaller aggregates for increasing φ_{σ}. For equal φ_{σ}, the slopes are shallower for the higher velocity. The powerlaw fit parameters can be found in Table 3. 
As shown in Fig. 9, the mass of the largest member of the powerlaw population, given by μ_{pw}, at first increases for larger values of φ_{σ}. It has a maximum at φ_{σ} ~ 0.02 for v_{0} = 10 m s^{1} and at φ_{σ} ~ 0.01 for v_{0} = 12.5 m s^{1}, and then decreases in an exponential fashion. When comparing both velocity cases (see Table 2), we find that more fragments are produced in the high velocity case. This was also observed by Geretshauser et al. (2011). The influence of the inhomogeneity on the powerlaw population can be summarised as follows: (1) the overall mass of the powerlaw distribution is increased as discussed above; (2) a higher degree of inhomogeneity leads to the production of smaller fragments; (3) the largest member of the powerlaw population reaches its highest mass at small φ_{σ}, after which its mass decreases with increasing φ_{σ}.
Fig. 8 Fragmentation parameter κ of the powerlaw population. For increasing inhomogeneity, indicated by the standard deviation φ_{σ}, the fragmentation parameter decreases. This indicates shallower slopes of the powerlaw mass distribution which is interpreted as an increasing fraction of smaller fragments. For φ_{σ} ≥ 0.02, the higher impact velocity results in smaller values of κ. 
Fig. 9 Fragmentation strength of the powerlaw population. The quantity μ_{pw} represents the normalised mass of the most massive member of the powerlaw population and serves as an indicator of the fragmentation strength. For increasing inhomogeneity, represented by larger values of φ_{σ}, μ_{pw} at first increases and then decreases again. The increase may be caused by our low quality statistics owing to the small number of fragments. 
Simulation results presented according to the fourpopulation classification.
Fit values of the powerlaw population.
5. Discussion and outlook
We have presented the first simulations of inhomogeneous porous preplanetesimals, where the inhomogeneity of SiO_{2} dust aggregates is explicitly resolved. This inhomogeneity can in principal also be determined in laboratory experiments (Güttler et al. 2009). The approach is based on the concept that according to our porosity model (Geretshauser et al. 2010) inhomogeneities in filling factor cause fluctuations in the compressive, shear, and tensile strength in the aggregate. These fluctuations can be regarded as flaws in the material. In contrast to damage models designed for brittle material (Grady & Kipp 1980), we did not explicitly evolve the propagation of these flaws. Instead, the defects in the dust material, which behaves more like a fluid, were approximated by the time evolution of the filling factor or – equivalently – the density. The inhomogeneity of an aggregate was assumed in the initial SPH particle distribution by using a Gaussian distribution to model the filling factor. The inhomogeneity is characterised by the standard deviation φ_{σ} of the Gaussian and the typical radius R_{c} of a clump.
Using this inhomogeneity approach, we have performed test simulations of collisions between dust aggregates of intermediate porosity with fixed typical clump size. Two collision velocities have been chosen: one below and one above the fragmentation threshold for homogeneous aggregates. The results have been analysed using the fourpopulation model of Geretshauser et al. (2011). For the lower collision velocity, we have shown that inhomogeneity leads to fragmentation. For both velocities, the masses of largest and second largest fragment are lower for a higher degree of inhomogeneity whereas the masses of the powerlaw and subresolution population are higher. Focussing on the powerlaw population, the number of fragments and the fraction of small fragments increase with increasing φ_{σ}. These findings show that the inhomogeneity approach has passed the first test in explicitly resolving the heterogeneity of porous objects. In a future study, we will address the influence of the typical clump size on the outcome of preplanetesimal collisions.
Our findings indicate that inhomogeneous dust aggregates are weaker than their homogeneous equivalents. A slight inhomogeneity is sufficient to cause a catastrophic disruption instead of growth as the result of a dust aggregate collision. Therefore, inhomogeneity might explain the higher velocity thresholds for fragmentation in simulations than the a lower value found in laboratory experiments (Geretshauser et al. 2011).
Furthermore, macroscopic dust aggregates in protoplanetary discs are produced by the subsequent impacts of smaller aggregates at different impact velocities, i.e. preplanetesimals have a collision history, thus are very likely to be inhomogeneous. This has a negative effect on planetesimal formation by subsequent collisions because preplanetesimals might become weaker as the number of collisions increases. As a result, their threshold velocity for catastrophic disruption might be smaller. However, in a future study we will investigate the influence of the clump size of the stability of the preplanetesimals. If the clumps are small relative to the preplanetesimal size, this might increase the stability of the aggregates.
With the model presented in this paper, we are able to include these features. Further studies might be carried out that investigate the inhomogeneity created by subsequent multiple impacts. The result could be classified according to the inhomogeneity model that we have presented here.
To date it has been found that the simulations of collisions between homogeneous aggregates carried out by means of the porosity model by Geretshauser et al. (2010) depend only on the filling factor and produce the same fragment distribution for an identical set of input parameters. By randomly assigning an inhomogeneity pattern it is now possible to profoundly investigate the statistics of a fragment distribution. Statistical fluctuations in the quantities of the fourpopulation model can now be estimated in particular for simulations with a small number of fragments.
In high velocity grazing collisions, the target begins to rotate. For highly porous aggregates, which have a low tensile strength, high spinning rates tend to lead to fragmentation of the target. With increasing inhomogeneity, this might also be true for aggregates with medium and low porosity. A quantitative investigation of this effect can be carried out by means of the presented approach.
The filling factor distribution of dust aggregates can be determined in the laboratory by Xray tomography measurements (Güttler et al. 2009). These empirical data can be directly implemented into our inhomogeneity model. In addition the successful material calibration (Geretshauser et al. 2010), the obtained inhomogeneity parameters can be used to improve the realistic simulation of porous dust aggregates.
Acknowledgments
R.J.G. and R.S. wish to thank S. Arena, W. Benz, and F. Meru for many illuminating discussions. We are obliged to our referee Gareth Collins for his profound revision and constructive comments which helped to improve and clarify the manuscript. The SPH simulations were performed on the university cluster of the computing centre of Tübingen, the bwGriD clusters in Karlsruhe, Stuttgart, and Tübingen. Computing time was also provided by the High Performance Computing Centre Stuttgart (HLRS) on the national supercomputer NEC Nehalem Cluster under project grant SPHPPC/12848. This project was funded by the Deutsche Forschungsgemeinschaft within the Forschergruppe “The Formation of Planets: The Critical First Growth Phase” under grant Kl 650/81.
References
 Benz, W. 2000, Space Sci. Rev., 92, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Benz, W., & Asphaug, E. 1994, Icarus, 107, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Benz, W., & Asphaug, E. 1995, Comp. Phys. Commun., 87, 253 [Google Scholar]
 Birnstiel, T., Ricci, L., Trotta, F., et al. 2010, A&A, 516, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 555 [Google Scholar]
 Geretshauser, R. J., Speith, R., Güttler, C., Krause, M., & Blum, J. 2010, A&A, 513, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Geretshauser, R., Meru, F., Speith, R., & Kley, W. 2011, A&A, 531, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 [Google Scholar]
 Grady, D. E., & Kipp, M. E. 1980, Int. J. Rock Mech. Mining Sci. Geomech. Abstr., 17, 147 [CrossRef] [Google Scholar]
 Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Herrmann, W. 1969, J. Appl. Phys., 40, 2490 [NASA ADS] [CrossRef] [Google Scholar]
 Hipp, M., & Rosenstiel, W. 2004, in EuroPar, ed. M. Danelutto, M. Vanneschi, & D. Laforenza (Springer), Lect. Notes Comp. Sci., 3149, 189 [Google Scholar]
 Jutzi, M., Benz, W., & Michel, P. 2008, Icarus, 198, 242 [NASA ADS] [CrossRef] [Google Scholar]
 Jutzi, M., Michel, P., Hiraoka, K., Nakamura, A. M., & Benz, W. 2009, Icarus, 201, 802 [NASA ADS] [CrossRef] [Google Scholar]
 Libersky, L. D., & Petschek, A. G. 1991, in Advances in the FreeLagrange method: including contributions on adaptive gridding and the smooth particle hydrodynamics method, ed. H. Trease, M. J. Fritts, & W. P. Crowley (Springer), Lect. Notes Phys., 395 [Google Scholar]
 Libersky, L. D., Petschek, A. G., Carney, T. C., Hipp, J. R., & Allahdadi, F. A. 1993, J. Comput. Phys., 109, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., & Stewart, G. R. 1993, in Protostars and planets, ed. E. H. Levy & J. I. Lunine (Univ. of Arizona Press), Space Sci. Ser., 1061 [Google Scholar]
 Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J. 2005, Rep. Progr. Phys., 68, 1703 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J., & Lattanzio. 1985, A&A, 149, 135 [NASA ADS] [Google Scholar]
 Randles, P. W., & Libersky, L. D. 1996, Comp. Meth. Appl. Mech. Engin., 139, 375 [Google Scholar]
 Schäfer, C. 2005, Ph.D. Thesis, Universität Tübingen, Tübingen [Google Scholar]
 Schäfer, C., Speith, R., & Kley, W. 2007, A&A, 470, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sirono, S. 2004, Icarus, 167, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Weibull, W. 1939, Ingeniörsvetenskapsakademiens handlingar, A statistical theory of the strength of materials (Stockholm: Generalstabens litografiska anstalts förlag), 151 [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and planets, ed. E. H. Levy, & J. I. Lunine (Univ. of Arizona Press), Space Sci. Ser., 1031 [Google Scholar]
 Wilkins, M. 1964, in Methods in Computational Physics V.3, ed. B. Alder, S. Fernbach, & M. Rotenberg (New York: Academic Press), 211 [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: List of symbols
List of symbols.
All Tables
Selection of the standard deviations φ_{σ} as they are used in the inhomogeneity algorithm.
All Figures
Fig. 1 Initial filling factor distributions of the target and projectile as binned (top) and cumulative diagram (bottom). The number density n(φ) of volumes with filling factor φ is plotted over φ. A larger degree of inhomogeneity is characterised by the increasing standard deviation φ_{σ} of the Gaussian (Eq. (17)) around the initial filling factor φ_{i} of the homogeneous aggregate. Here, we choose φ_{i} = 0.35. The values of φ_{σi} can be found in Table 1. 

In the text 
Fig. 2 Interior and exterior view of the targets for different inhomogeneities with the filling factor colour coded. The inhomogeneity increases from the homogeneous target a) to the most inhomogeneous target f). In particular, the standard deviations are φ_{σ} = 0 a), φ_{σ} = 0.01 b), φ_{σ} = 0.02 c), φ_{σ} = 0.03 d), φ_{σ} = 0.04 e), and φ_{σ} = 0.05 f). All targets have a similar pattern that originates from considering interacting particle packages in the implementation. With increasing φ_{σ} the maximum and minimum filling factor of the different spots increase and decrease, respectively. The projectiles have a similar appearance. The result of collisions among these aggregates with 10 m s^{1} are depicted in Fig. 3. 

In the text 
Fig. 3 Outcome of a collision between a target and projectile with the same φ_{σ} for different inhomogeneities. In all cases, the target and projectile radii are r_{t} = 10 cm and r_{p} = 6 cm, respectively, and the collision velocity is v_{0} = 10 m s^{1}. Analogous to the initial targets in Fig. 2, the standard deviations for the Gaussian of the initial targets are φ_{σ} = 0 a), φ_{σ} = 0.01 b), φ_{σ} = 0.02 c), φ_{σ} = 0.03 d), φ_{σ} = 0.04 e), and φ_{σ} = 0.05 f). The collision outcome is shown in the impact direction. In the homogeneous case a) and for small inhomogeneities b), the target stays intact and forms one massive object with the projectile. For φ_{σ} ≥ 0.02, the target fragments c)–f). The fragment sizes decrease with increasing φ_{σ} and at the same time the number of fragments increases. 

In the text 
Fig. 4 The total masses of the four populations with inhomogeneities measured by φ_{σ}. The collision velocities are v_{0} = 10 m s^{1} (top) and v_{0} = 12.5 m s^{1} (bottom). The masses of the largest fragment m_{1}, of the second largest fragment m_{2}, of the powerlaw population m_{pw}, and of the subresolution population m_{sr} are stacked and normalised by the total mass of the system m_{tot}. In both cases, m_{1} and m_{2} decrease with increasing φ_{σ}. Conversely, m_{pw} and m_{sr} increase. For the homogeneous case (φ_{σ} = 0), sticking is found for the low velocity collision, whereas for v_{0} = 12.5 m s^{1}m_{2} ≠ 0 and m_{pw} ≠ 0, which indicates fragmentation. 

In the text 
Fig. 5 The energy of each population for increasing inhomogeneity parameter φ_{σ} and collision velocity v_{0} = 10 m s^{1} (top) and v_{0} = 12.5 m s^{1} (bottom). The energy of each population is the sum of translational, rotational and vibrational energy and normalised by the initial kinetic energy E_{init}. The two latter energy contributions are negligible for the headon collisions presented here. Owing to the strong correlation with mass, this diagram follows a similar trend to Fig. 4. The energy contributions are divided up into E_{1}, E_{2}, E_{pw}, and E_{sr} for the largest and second largest fragment, the powerlaw and subresolution population, respectively. In both velocity cases, both E_{1} and E_{2} decrease with increasing φ_{σ}. In contrast, E_{pw} and E_{sr} increase. In particular, the high fraction of kinetic energy stored in the subresolution population is remarkable. 

In the text 
Fig. 6 Cumulative mass of the fragment masses over their average filling factor φ. The mass is normalised by the total fragment mass m_{tot}. The collision velocities are v_{0} = 10 m s^{1} (top) and v_{0} = 12.4 m s^{1} (bottom). Compared to the initial distribution of filling factors (Fig. 1), the curves are shifted significantly towards higher filling factors. This reflects the compression taking place during the collision. 

In the text 
Fig. 7 Cumulative mass distribution of a powerlaw population for different inhomogeneity parameters φ_{σ} and collision velocities v_{0} = 10 m s^{1} (top) v_{0} = 12.5 m s^{1} (bottom). In both velocity cases both the mass of the largest member of the powerlaw population and the powerlaw slope decrease. This indicates the fragmentation to smaller aggregates for increasing φ_{σ}. For equal φ_{σ}, the slopes are shallower for the higher velocity. The powerlaw fit parameters can be found in Table 3. 

In the text 
Fig. 8 Fragmentation parameter κ of the powerlaw population. For increasing inhomogeneity, indicated by the standard deviation φ_{σ}, the fragmentation parameter decreases. This indicates shallower slopes of the powerlaw mass distribution which is interpreted as an increasing fraction of smaller fragments. For φ_{σ} ≥ 0.02, the higher impact velocity results in smaller values of κ. 

In the text 
Fig. 9 Fragmentation strength of the powerlaw population. The quantity μ_{pw} represents the normalised mass of the most massive member of the powerlaw population and serves as an indicator of the fragmentation strength. For increasing inhomogeneity, represented by larger values of φ_{σ}, μ_{pw} at first increases and then decreases again. The increase may be caused by our low quality statistics owing to the small number of fragments. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.