Issue 
A&A
Volume 587, March 2016



Article Number  A128  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201526565  
Published online  02 March 2016 
Comet formation in collapsing pebble clouds
What cometary bulk density implies for the cloud mass and dusttoice ratio
^{1}
MaxPlanck Institute for Solar System Research,
JustusvonLiebig Weg 3,
37077
Göttingen,
Germany
email:
lorek@mps.mpg.de
^{2}
Institut für Geophysik und extraterrestrische Physik, Technische
Universität Braunschweig, Mendelssohnstr. 3, 38106
Braunschweig,
Germany
Received: 20 May 2015
Accepted: 13 January 2016
Context. Comets are remnants of the icy planetesimals that formed beyond the ice line in the solar nebula. Growing from μmsized dust and ice particles to kmsized objects is, however, difficult because of growth barriers and time scale constraints. The gravitational collapse of pebble clouds that formed through the streaming instability may provide a suitable mechanism for comet formation.
Aims. We study the collisional compression of silica, ice, and silica/icemixed pebbles during gravitational collapse of pebble clouds. Using the initial volumefilling factor and the dusttoice ratio of the pebbles as free parameters, we constrain the dusttoice mass ratio of the formed comet and the resulting volumefilling factor of the pebbles, depending on the cloud mass.
Methods. We use the representative particle approach, which is a Monte Carlo method, to follow cloud collapse and collisional evolution of an ensemble of ice, silica, and silica/icemixed pebbles. Therefore, we developed a collision model which takes the various collision properties of dust and ice into account. We study pebbles with a compact size of 1 cm and vary the initial volumefilling factors, φ_{0}, ranging from 0.001 to 0.4. We consider mixed pebbles as having dusttoice ratios between 0.5 and 10. We investigate four typical cloud masses, M, between 2.6 × 10^{14} (very low) and 2.6 × 10^{23} g (high).
Results. Except for the very lowmass cloud (M = 2.6 × 10^{14} g), silica pebbles are always compressed during the collapse and attain volumefilling factors in the range from ⟨ φ ⟩ _{V} ≈ 0.22 to 0.43, regardless of φ_{0}. Ice pebbles experience no significant compression in very lowmass clouds. They are compressed to values in the range ⟨ φ ⟩ _{V} ≈ 0.11 to 0.17 in low and intermediatemass clouds (M = 2.6 × 10^{17}−2.6 × 10^{20} g); in highmass clouds (M = 2.6 × 10^{23} g), ice pebbles end up with ⟨ φ ⟩ _{V} ≈ 0.23. Mixed pebbles obtain filling factors in between the values for pure ice and pure silica. We find that the observed cometary density of ~0.5 g cm^{3} can only be explained by either intermediate or highmass clouds, regardless of φ_{0}, and also by either very low or lowmass clouds for initially compact pebbles. In any case, the dusttoice ratio must be in the range of between 3 ≲ ξ ≲ 9 to match the observed bulk properties of comet nuclei.
Key words: comets: general / methods: numerical
© ESO, 2016
1. Introduction
Comets are the remnants of icy planetesimals that formed beyond the snow line in the outer part of the solar nebula about 4.6 Gyr ago. They typically have sizes from 1 km to 10 km (Lamy et al. 2004) and consist of different ices (mainly water ice) and dust. The mean density of comet nuclei ranges from 0.4 ± 0.2 g cm^{3} (A’Hearn 2011) to 0.6 ± 0.4 g cm^{3} (Blum et al. 2006). Recently, Sierks et al. (2015) determined the mean density of comet 67P/ChuryumovGerasimenko at 0.470±0.045 g cm^{3}. Comparing these estimates, ρ_{c} = 0.5 ± 0.3 g cm^{3} can be considered as a typical comet density. Therefore, to obtain such a low density, comets must be porous objects with a nucleus structure that can be thought of as a rubblepilelike assemblage of dust and ice, held together by selfgravity rather than material strength (Weissman et al. 2004).
Growth from μmsized dust and ice particles to kmsized comets is difficult to explain because the process is affected by several growth barriers: the bouncing barrier for mm to cmsized particles (Zsom et al. 2010), the fragmentation barrier for cm to msized bodies (Blum & Wurm 2008; Güttler et al. 2010), and the radial drift barrier, which peaks for msized objects at a heliocentric distance of ~1 au (Weidenschilling 1977).
One possible scenario for overcoming the growth barriers employs the bouncing barrier. A big particle sweeps up the small particles stuck at the bouncing barrier and thereby grows in masstransfer collisions (Windmark et al. 2012a; Dra¸żkowska et al. 2013). It was found that a 100 msized body can form at a heliocentric distance of 3 au within 1 Myr. However, the presence of big seed particles is necessary to trigger this mechanism. Among the main sources of collision velocities of small dust grains in the protoplanetary disk are Brownian motion and turbulence, which are stochastic processes. Thus, it is reasonable to assume that the collision velocities are distributed around a typical collision speed. This velocity distribution then helps to form those seeds, because some particles fortuitously undergo only lowvelocity sticking collisions and thus grow past the barriers (Windmark et al. 2012b,c). However, the dating of calciumaluminiumrich inclusions (CAIs) and chondrules found in meteorites puts constraints on the timescale for the formation of planetesimals inside the snow line, which turns out to be a time span of only a few Myr after solar system formation (Connelly et al. 2012). For the cometesimals outside the snow line, the formation period is less constrained, but the lifetime of the gaseous disk of ≲10 Myr (Fedele et al. 2010) is an upper limit. Therefore, the formation of kmsized bodies via sweepup growth is not fast enough, especially at large semimajor axes, where growth is generally slower than in the inner solar system.
The streaming instability (SI) offers another mechanism for planetesimal formation. This type of twofluid instability arises in the protoplanetary disk because dust and gas rotate at different velocities as a result of the subKeplerian speed of the gas and the nonperfect coupling of the pebbles to the gas in the disk. Using a linear perturbation analysis of the hydrodynamic equations, Youdin & Goodman (2005) showed that this situation is intrinsically unstable and tends to cluster dust particles, eventually leading to the formation of planetesimals.
The SI avoids the growth barriers by jumping directly from pebbles with Stokes number St ~ 0.01−1 (corresponding to a grain size of mm to dm, depending on the radial location in the disk) to Ceressized planetesimals (Johansen et al. 2007; Bai & Stone 2010; Dra¸żkowska & Dullemond 2014). The SI is triggered once the local dusttogas ratio in the protoplanetary disk exceeds unity. In this case, the feedback of the dust on the gas motion via the drag force causes the local dust density enhancement to grow in a runaway process. This way, the local dust density grows up to 1000 times the local gas density. If the density reaches Roche density, a gravitationally bound clump is formed. Subsequent gravitational collapse can finally form a planetesimal (Johansen et al. 2007, 2014; Wahlberg Jansson & Johansen 2014).
As argued by Skorov & Blum (2012) and Blum et al. (2014), the benefits of SI as a mechanism for forming comets is that it bypasses growth barriers, while the collision speeds of pebbles in a collapsing cloud also remain low, at most the escape velocity of the cloud. Hence, pebble fragmentation plays only a minor role, if at all, and the dominant collisiontype is bouncing. The collapse of a cloud of bouncing pebbles then leads to a comet with a packing fraction of pebbles of φ_{p} ≈ 0.6. This packing fraction means that, assuming the comet is a sphere composed of spherical pebbles, only 60% of the comet’s volume is actually filled with pebbles; the remaining 40% is void space. However, a packing fraction of 0.6 assumes equalsized pebbles. A size distribution changes the packing fraction, because the voids created by the big pebbles can be be filled with small pebbles, which results in a denser packing. The pebbles themselves are also not compact spheres, but possess a filling factor of φ ≈ 0.4. As with the packing fraction, the volumefilling factor denotes that only 40% of the pebble’s volume is filled with solid material (ice or dust); again, the remaining 60% is void space. This fluffiness is a result of the growth process in the protoplanetary disk. Hitandstick collisions of μmsized dust grains form fractal aggregates with lowfilling factors. After reaching a certain size of typically a few mm, the dust aggregates encounter the bouncing barrier and get compressed in bouncing collisions, which leads to a filling factor of ~0.4 (Zsom et al. 2010). Experimental work by Weidling et al. (2009) and Güttler et al. (2009) also showed that the volumefilling factor of silica pebbles acquires a value of ~0.4 after experiencing many bouncing collisions. The filling factor of the comet, φ_{c}, is then the product of both: φ_{c} = φ_{p} × φ ≈ 0.24 (Skorov & Blum 2012). Comets are thus expected to be porous bodies, which agrees with recent density measurements and estimates (e.g. Blum et al. 2006; A’Hearn 2011; Sierks et al. 2015).
As previously mentioned, ice is a significant constituent of comets. Collision models based on laboratory experiments (Güttler et al. 2010; Windmark et al. 2012a) so far only deal with silica aggregates, because experiments with ice particles are difficult to conduct. However, early results suggest that the threshold velocities for sticking, bouncing, and fragmentation of ice particles are by a factor of 10 higher than for silica (Gundlach et al. 2011; Gundlach & Blum 2015; Hill et al. 2015). The compression curve of ice aggregates, on the other hand, cannot be scaled as easily as the threshold velocities (see Fig. 1).
Here we describe numerical simulations of collapsing pebble clouds formed through the SI that are based on the energy formalism developed by Wahlberg Jansson & Johansen (2014). We implement an improved collision model which allows us to use dust, ice, and dust/icemixed pebbles and to follow the filling factor of these pebbles. We conduct three sets of simulations, one with silica, one with ice, and one with mixed pebbles, to derive the compression behaviour of the pebbles during the collapse. Furthermore, we use the observed density of cometary nuclei and the estimated pebble packing to constrain the dusttoice ratio of the mixed pebbles, the cloud mass, and the pebble properties.
The paper is structured as follows: we briefly describe the laboratory experiments that lead to the compression curves of ice aggregates in Sect. 2. Then, we outline the numerical method in Sect. 3. We present the simulation outcomes in Sect. 4. In Sect. 5, we discuss the results and the uncertainties. Finally, we conclude the paper by summarising our findings in Sect. 6.
2. Laboratory experiments on dust and ice compression
Compression experiments on dust samples have been performed by Blum & Schräpler (2004); Blum et al. (2006); Güttler et al. (2009) using silica particles as dust analogues. Güttler et al. (2009) approximate the compression curves by an analytic function of the form (1)Here, p is the applied pressure, p_{m} is the pressure at which 50% of the compression is reached, and Δ is the width of the transition regime. The factors, φ_{1} and φ_{2} are the minimum and the maximum volumefilling factors of the material, respectively. Depending on the nature of the compression experiments, i.e. unidirectional or omnidirectional, the maximum volumefilling factor varies between φ ≈ 0.3 (unidirectional) and φ ≈ 0.6 (omnidirectional). Table 1 summarises the results of their experiments on dust compression.
Güttler et al. (2009) give a value of p_{m} = 1.3 × 10^{5} dyn cm^{2} for the static omnidirectional compression of the dust samples. To fit their dynamic compression experiments, they refined this value to 1.3 × 10^{4} dyn cm^{2} by comparing their measurements to the results of smoothedparticle hydrodynamics (SPH) simulations. Güttler et al. (2010) adopt the lower value in their model for the porosity evolution and furthermore extrapolate Eq. (1)for low pressures (p<p_{m}) with a powerlaw function, (2)to account for volumefilling factors lower than the ones used in their experiments. Hence, we stick with the lower value of p_{m} and a compression curve for silica given by Eqs. (1)and (2)(see dotted line in Fig. 1).
Fit parameters for the silica and the granular water ice samples.
In contrast to dust, the compression properties of granular water ice have, so far, not been measured. Thus, we performed laboratory experiments to investigate how the volumefilling factor of granular water ice changes when an external load is applied to the sample.
We produce granular water ice samples by spraying liquid water into a cooled sedimentation chamber (this method has been extensively discussed in previous works, see e.g. Gundlach et al. 2011; Jost et al. 2013). This method produces spherical micrometresized water ice particles with a mean radius of 1.45 μm. After sedimentation of the micrometresized water ice particles in the sedimentation chamber, a porous water ice aggregate with an initial volumefilling factor of ~0.1 is formed on a cold plate at the bottom of the chamber.
Fig. 1 Compression curves of the granular water ice samples consisting of micrometresized water ice particles (coloured curves). The initial volumefilling factor of the samples varies between 0.08 and 0.12. The dashed curve shows the best fit to the ice measurements. For comparison, the dotted line is the omnidirectional compression curve for silica particles used in the simulations (Güttler et al. 2009, 2010). 
These samples are then transported to a second experimental setup that was constructed to measure the compressive strength of the water ice samples. Here, the samples are placed, together with the cold plate, inside a transport container filled with vaporised nitrogen to avoid sintering^{1} and condensation of frost on the sample’s surface. The typical temperature of the cold plate during the transport process is ~100 K and the transfer time is less than 30 s.
To perform the compressive strength measurements, the transport container is placed on top of a scale and beneath a cooled piston (the temperature of the cold plate and the piston have never exceeded 125 K during the experiments). Then, the piston is slowly pressed into the granular sample by a stepper motor with a speed of 10 μm s^{1}. The piston exerts a local compression on the sample, i.e. the material is able to flow into the noncompressed volumes of the sample.
A total of four measurements are performed with this experimental setup (solid curves in Fig. 1). The initial volumefilling factor varies between 0.08 and 0.12. Relatively low compressions (<10^{3} dyn cm^{2}) do not lead to an increase in the volumefilling factor of the water ice samples. However, higher compressions are sufficient to compress the granular water ice samples from the initial volumefilling factor to a volumefilling factor of ~0.27 (for a compression of 10^{6} dyn cm^{2}). The experiments stop at 10^{6} dyn cm^{2}, because the existing experimental setup (we adopted the experimental setup that was used for the dust measurements by Güttler et al. 2009) is not capable of detecting higher pressure values.
For comparison, the compression curve of silica dust aggregates composed of micrometresized silica particles is shown by the dotted curve in Fig. 1 (see Güttler et al. 2009, 2010, for details). The comparison of the compression curves directly shows that the water ice aggregates require more load than the dust aggregates to reach the same volumefilling factor. This effect can be explained by the higher specific surface energy of water ice, compared to the silica dust particles (see Gundlach et al. 2011). Compression of granular matter in this pressure regime is caused by the relocation of the particles in the sample by the rolling of the particles into the voids of the material, which then induces an increase of the mean coordination number per particle. Because of the higher specific surface energy of water ice, the rolling friction force of water ice is also enhanced in comparison to the silica dust particles, which implies that a higher pressure is required to rearrange the water ice particles in the samples. For a more detailed discussion of the compression physics of granular material, see Schräpler et al. (2015).
Since the granular ice aggregates are harder to compress, the experimental setup is not capable of resolving the entire pressure curve. Thus, the pressure curves of the water ice samples do not saturate, as it is the case for comparable dust aggregates. Filling factors higher than 0.27 might thus be possible for pressures exceeding 10^{6} dyn cm^{2}, but more experiments with a modified setup would be required to test this and, thus, outside the objectives of this paper.
To analyse the data in more detail, we apply Eq. (1)based on the Fermi distribution with logarithmic pressure to fit each individual data set. The values of the fit parameters are summarised in Table 1. The application of the Fermi distribution is motivated by the fact that Güttler et al. (2009) use a similar function to fit the compression curves of the silica dust aggregates. To provide a general compression curve for granular water, we calculated the geometric mean of the individual fit parameters (see Table 1).
A comparison of the fit parameters of the water ice samples with the silica dust aggregates (unidirectional and omnidirectional, Table 1) indicates that the water ice samples behave similarly to the unidirectional compression of the dust aggregates (see Güttler et al. 2009, for details). This means that the material is able to flow into the noncompressed areas of the samples. Only the parameter Δ is too high for fitting into this picture. This indicates that the compressional behaviour of water ice cannot only be explained by shifting the compression curves of the dust aggregates to higher compression pressures.
3. Numerical method
A comet of radius 5 km forming from ice pebbles with 1 cm radius, for example, would require ~10^{17} pebbles – a number that is impossible to simulate directly. Therefore, we employ the representative particle approach, developed by Zsom & Dullemond (2008), for dust coagulation in the protoplanetary disk to model the collapse of the pebble cloud. The Monte Carlo approach is suitable because it allows modelling a large number, N, of physical particles by using only a small random sample of those particles. This means, n ≪ N particles out of the N particles are randomly chosen, which are the representative particles. Thus, the fundamental assumption is that the distribution function of the representative particles can be regarded as the distribution function of the physical particles. The time evolution of the whole system is then given by following the change of the properties of the n representative particles over time. The random sample can be small, which has the advantage of being computationally feasible, e.g. 100 representative particles in contrast to the 10^{17} physical particles.
3.1. Cloud collapse model
The details of the collapse model are outlined in Wahlberg Jansson & Johansen (2014). The model employs energy considerations for a pebble cloud that starts in virial equilibrium. Without any means of energy dissipation, the cloud remains in equilibrium because the pebbles move on orbits oscillating between apo and pericentre. Elastic collisions only redistribute the orbital planes, which leads to random motion (pressure) counteracting gravity. The collapse sets in, because energy is dissipated in inelastic pebble collisions. The cloud contracts and releases gravitational energy to compensate for the dissipated energy. The collapse ends, once the cloud reaches the desired density at which all pebbles clump together and a solid body is formed. In our simulations, the collapse stops when the cloud reaches a density of 0.5 g cm^{3}, which we consider here as the typical density of a comet.
3.2. Collision model for dust and ice pebbles
To simulate collisions between pebbles, we use the collision model by Windmark et al. (2012a). Based on the extensive work by Güttler et al. (2010), Windmark et al. (2012a) identify four collision types as the most important for dust growth in the protoplanetary disk: sticking, bouncing, fragmentation, and partial fragmentation with mass transfer. This collision model is supported by laboratory experiments of colliding dust aggregates and provides mass and velocity dependent thresholds on the four collision regimes. Hence, it can be applied not only to modelling the dust growth in the protoplanetary disk, but also to model the pebble collisions in the collapsing cloud.
The threshold velocities for sticking, bouncing, and fragmentation of dust aggregates are given by Windmark et al. (2012a): where m_{p,t} is the mass of the projectile, or target, respectively. The normalising constants m_{s} = 3.0 × 10^{12} g and m_{b} = 3.3 × 10^{3} g are determined using laboratory data. The normalisation mass m_{1.0} = 3.67 × 10^{7} g describes the onset of fragmentation. By replacing m_{1.0} with m_{0.5} = 9.49 × 10^{11} g, one obtains the velocity threshold for which the largest remnant is half of the mass of the projectile or target. The fragmentation threshold velocity is calculated in the centre of mass reference frame. Thus, Δv_{frag} is compared to the centre of mass velocities of target and projectile, respectively: For equalmass pebbles, the centre of mass velocity is half of the relative velocity.
Comets are very porous objects. Hence, we add the porosity evolution recipe for sticking and bouncing described in Güttler et al. (2010) to the collision model. Sticking collisions form fractal aggregates and the volumefilling factor decreases. Ormel et al. (2007) describe how to calculate the new filling factor of the resulting aggregate. In contrast, bouncing collisions compress the pebbles. The maximum compression per collision depends on collision velocity and pebble density, and thus on the pressure exerted on target and projectile. The compression curve gives the maximum filling factor that can be obtained for a certain pressure (see Fig. 1). The details are outlined in Güttler et al. (2009, 2010) and Weidling et al. (2009). However, for fragmentation and mass transfer, we simplistically assume no porosity change.
To include water ice pebbles into the simulations, we scale the threshold velocities Eqs. (3)−(5)by a factor of ten to higher values, i.e. , and similarly for bouncing and fragmentation. This factor of ten is motivated by the fact that the specific surface energy and the sticking threshold velocity of water ice is also ten times higher (Gundlach et al. 2011; Gundlach & Blum 2015; Hill et al. 2015). The measured compression curve of ice shown in Fig. 1 does not allow for an easy scaling. Thus, we use Eq. (1)with the mean values of the fit parameters given in Table 1.
3.3. Collision model for dust/icemixed pebbles
Using pure ice and pure silica pebbles is a major simplification. There is evidence for μmsized pure water ice grains without refractory inclusions in the interior of comets (Sunshine et al. 2007). However, comet nuclei are also known to contain crystalline silicates, which cannot have formed in the cold outer part of the solar nebula (McKeegan et al. 2006). Thus, radial mixing must have taken place in the solar nebula (BockeléeMorvan et al. 2002; Hanner & Bradley 2004; Cuzzi & Weidenschilling 2006). The coagulation timescale to grow from μmsized grains to cmsized pebbles is expected to be short (≲10^{3} yr, Zsom et al. 2010) compared to the radial mixing timescale (≳10^{4} yr, BockeléeMorvan et al. 2002).
Therefore, we assume that small silica and ice grains (monomers) originate inside and outside the snow line, respectively, and are mixed together in the cometforming region. Thereby, the monomers coagulate to form cmsized dust/icemixed pebbles. Unfortunately, we do not know how mixed pebbles behave in collisions. The outcome will likely be in between the pure silica and the pure ice case, depending on the dust/ice content. We propose a simple interpolation scheme between both cases for mixed pebbles, which is based on the fractional abundances of dust and ice monomers (see Sect. 3.3.2).
For our interpolation scheme to work, we assume a homogeneous mixture of silica and ice monomers within the mixed pebble. A pebble with a refractory core and an icy mantle is expected to behave like an icy pebble, but with the mass of a silica pebble, because the density ratio of silica to ice is ~3. However, we do not consider coremantle pebbles in our study.
3.3.1. Density of a mixed pebble
A single pebble of mass m is composed of N_{d} dust monomers of mass m_{0,d} and N_{i} ice monomers of mass m_{0,i}. The total mass of dust and ice within the pebble is then m_{d} = N_{d}m_{0,d} and m_{i} = N_{i}m_{0,i}, respectively, and the pebble composition is expressed in terms of the dusttoice ratio, ξ ≡ m_{d}/m_{i}. The dusttoice ratio is the mass ratio of dust to ice and not the volume ratio. Furthermore, we only consider water ice and neglect other ice species, such as carbon monoxide, carbon dioxide, methane, or ammonia.
We first calculate the density of a mixed pebble with fixed dusttoice ratio ξ. The pebble mass is the sum of the masses of all dust and ice monomers: (8)where V_{d,i} are the compact volumes occupied by dust and ice monomers, respectively, and ρ_{d,i} are the corresponding monomer densities. On the other hand, the pebble mass is given by (9)where φ is the pebblefilling factor and ρ_{•} is the pebble density. Using V = V_{d} + V_{i} and ν_{d,i} = V_{d,i}/V, it follows that ν_{d} + ν_{i} = 1. Together with the definition of ξ, Eqs. (8)and (9)can be solved for ρ_{•}: (10)where is the compact density without any void space between the monomers. Thus, the density of the pebble depends on the dusttoice ratio. For ξ = 0 (no dust), and for ξ = ∞ (no ice) we get , as expected.
3.3.2. Interpolation scheme
To interpolate the collision behaviour between the two cases of pure dust and pure ice, we use the fractional abundance of dust monomers, x = N_{d}/ (N_{d} + N_{i}), because the thresholds for sticking, bouncing, and fragmentation depend on the ability of the collision to rearrange or break bonds between the monomers (Dominik & Tielens 1997; Blum & Wurm 2000). If there are more dust monomers than ice monomers, the pebble behaves more like a dust pebble. In the opposite case, when there are more ice monomers than dust monomers, the pebble behaves more like ice.
Thus, we formulate the sticking threshold velocity for a mixed pebble as (11)Here, we use the fact that the threshold velocity for ice is ten times the value for dust (see Sect. 3.2). Similar expressions can be derived for the bouncing and fragmentation threshold velocity.
However, the remaining recipe of Windmark et al. (2012a) for determining masses of the remnant and fragments does not change in this case. However, there is one important assumption in the case of fragmentation: the composition of the pebble does not change. If, prior to the collision, the pebble has a dusttoice ratio of ξ, the remnant and the fragments will have the same value for the dusttoice ratio after the collision. This should not cause any problems in fragmentation/coagulation cycles, unless pebbles are ground down to monomers. In this case, the composition can change as a result of the accretion of single monomers. The collision velocities in the collapsing pebble cloud, however, are too low for this effect to play an important role.
Pebbles are compressed in bouncing collisions. The maximum compression is given by the compression curve φ(p) (see Eqs. (1)and (2)), which provides the maximum volumefilling factor, φ, that a pebble obtains for a given pressure p. As with the threshold velocities for sticking, bouncing, and fragmentation, the compression of a pebble is determined by rearranging bonds between monomers. Thus, we use the same linear interpolation scheme based on the fractional abundance of dust monomers, x, to get (12)where φ^{d,i}(p) are the compression curves of pure silica and pure ice aggregates, respectively, which cannot be simply scaled as is the case for threshold velocities (see Fig. 1 and Table 1).
The change of the volumefilling factor in sticking collisions (see Ormel et al. 2007, for details) is valid also for mixed pebbles. To compensate the suppression of fractal growth in collisions between aggregates and single particles, Ormel et al. (2007) introduce an additional term, Ψ_{add} ∝ exp(−μ/m_{F}) (μ is the reduced mass of target and projectile), to their Eq. (13). The mass scale m_{F} = 10 m_{0}, which ensures the correction term vanishes for collisions between two massive aggregates, depends on the monomer mass m_{0}. For mixed pebbles, however, silica and ice monomers are present and thus we replace m_{0} by the reduced monomer mass in our implementation of the Ormel et al. (2007) porosity model.
3.4. Initial conditions
With this setup (collapse and collision model), we conduct three sets of simulations. The initial compact pebble size, without void space, is s^{∗} = 1 cm. We vary the initial filling factor, φ_{0}, of the pebbles and the pebble cloud mass, M. In the first set, we use water ice pebbles with a compact density of . In the second set, we have silica pebbles with a compact density of . In the third set, we simulate mixed pebbles with varying dusttoice ratios, ξ, and densities according to Eq. (10).
The cloud mass equals the mass of the final cometesimal, because it is assumed that the cloud collapses into one single object (Wahlberg Jansson & Johansen 2014). We use four different cloud masses equivalent to the mass of an object with radius 0.5 km, 5 km, 50 km, and 500 km, respectively, with a bulk density of ρ_{c} = 0.5 g cm^{3}, which is formed by the collapse of the cloud. We refer to these four masses as very low, low, intermediate, and highmass cloud or case, respectively. The exact masses of the clouds are 2.6 × 10^{14} g (verylow mass), 2.6 × 10^{17} g (low mass), 2.6 × 10^{20} g (intermediate mass), and 2.6 × 10^{23} g (high mass). The smallest size represents the size of a comet, whereas the largest one corresponds to the size of Ceres.
We use initial filling factors in the range between 0.001−0.4 to span the gap from very fluffy aggregates to compact pebbles that form when reaching the bouncing barrier. We do not expect filling factors much higher than 0.4 because this is the limit given by the compression properties of the aggregates at the bouncing barrier (Weidling et al. 2009; Zsom et al. 2010). Very fluffy aggregates can form in hitandstick collisions of icy particles outside the snow line (e.g. Okuzumi et al. 2012; Kataoka et al. 2013). Filling factors lower than 0.001 are, in principle, possible in hitandstick collisions of ice (Okuzumi et al. 2012, including collisional compression of the combined aggregate) but, due to the predominantly bouncing, we can always expect a certain amount of compression, which increases the filling factor rapidly. Furthermore, in the simulations with mixed pebbles, we use four different initial filling factors ranging from 0.1 to 0.4 and six different dusttoice ratios between 0.5 (icedominated) and 10 (dustdominated). We use a constant coefficient of restitution of ≈0.7 (Blum & Münch 1993; Weidling et al. 2012; Weidling & Blum 2015). However, changing the value of the coefficient of restitution slightly affects the final filling factor of the pebbles. We address this issue in the caveats.
3.5. Resolution study
Fig. 2 Resolution study. Volumeweighted mean filling factor of pebbles as a function of the number of representative particles for different initial filling factors (φ_{0} = 0.01, 0.20, and 0.40) for the lowmass cloud. Increasing the number of particles does not change the final filling factor. 
We tested different numbers of representative particles and found no significant resolution dependency for the final volumefilling factor of the pebbles. Figure 2 shows the volumeweighted mean filling factor of ice and silica pebbles for an increasing number of representative particles in the lowmass case. Increasing the number of particles does not change the final filling factor because each representative particle experiences, on average, ≳2000 collisions. This is in accordance with Weidling et al. (2009) who reported saturation of the volumefilling factor after ~2000 collisions. The mean number of collisions per representative particle also exceeds 2000 in the intermediate and highmass cloud for n = 100, which leads, together with the higher collision velocities in these cases, to a volumefilling factor approaching the maximum value (0.46 for silica and 0.27 for ice according to Güttler et al. 2010, and Fig. 3). Thus, we do not expect any changes when increasing the number of particles. However, there is some scatter in the final collapse time, which can be reduced by using more particles but, on average, we successfully reproduce the results found by Wahlberg Jansson & Johansen (2014). Based on these tests, we use n = 100 representative particles for the runs presented here.
4. Results
In this section, we present the results of our two sets of simulation runs. In Table 2, we provide the values for the final filling factor of the pebbles obtained in our simulations as an overview.
4.1. Collision regimes
We find that bouncing is the most frequent collision type in all four cloud masses, accounting for almost 100% of pebble collisions. Some sticking occurs for ice, dust, and mixed pebbles alike in the very lowmass cloud and for small fragments in higher mass clouds. Fragmentation and mass transfer contribute to the collision outcomes for higher cloud masses, which leads to a shift in the size distribution towards smaller pebbles in these clouds. However, the contribution by sticking, fragmentation, or mass transfer is insignificant compared to the very large number of bouncing collisions. Thus, for studying the porosity evolution of the pebbles, the a priori assumption of neglecting the porosity changes in fragmenting collisions in our collision model is justified.
4.2. Volumefilling factor
Fig. 3 Pebble compression for pure ice and silica pebbles. Volumeweighted mean filling factor at the end of collapse as a function of initial filling factor: a) silica pebbles; b) ice pebbles. The four cases of a very lowmass (blue), a lowmass (yellow), an intermediatemass (green), and highmass (red) cometesimal are shown. Along the black dotted line, the final filling factor equals the initial filling factor. The error bar indicate the standard deviation of the pebble ensemble. Silica pebbles are significantly compressed for M> 2.6 × 10^{20} g. The maximum filling factor ice pebbles can obtain by compression is considerably lower than for silica, φ ~ 0.23 compared to φ ~ 0.43. 
Volumeweighted mean filling factor of pebbles for all runs performed in this study.
Fig. 4 Pebble compression for mixed pebbles. Volumeweighted mean filling factor at the end of collapse as a function of initial filling factor: a) very lowmass cloud; b) lowmass cloud; c) intermediatemass cloud; d) highmass cloud. Mixed pebbles with different dusttoice ratios are shown: 0.5 (blue, icedominated), 1.0 (yellow), 3.0 (green, equal amount of dust and ice), 5.0 (red), 7.0 (violet), and 10 (brown, dustdominated). Along the black solid line, the final filling factor equals the initial filling factor. The error bar indicate the standard deviation of the pebble ensemble. As for pure silica or ice pebbles, compression increases with increasing cloud mass. Furthermore, compression increases with increasing dusttoice ratio. Pebbles with low dusttoice ratios behave like ice pebbles, while pebbles with high dusttoice ratios behave more like silica pebbles. 
We investigate the final volumeweighted mean filling factor of pebbles: (13)where N_{i} = M/ (nm_{i}) is the number of particles a representative particle represents, is the volume of a pebble ( is the compact volume, i.e. the volume without void space), and φ_{i} is the volumefilling factor of a pebble. We weigh the filling factor with the volume, because we expect large and fluffy pebbles to dominate the porosity of the comet. The mean value in Eq. (13)is a weighted sample mean of the pebble ensemble. To quantify the uncertainty of the volumefilling factor, we use the weighted sample standard deviation of the pebbles.
4.2.1. Pure dust and ice pebbles
We find a clear difference in ⟨ φ ⟩ _{V} for ice and silica and for different cloud masses. As can be seen in Fig. 3a for the very lowmass cloud, compression of silica pebbles is negligible and the initial filling factor is retained, unless rare sticking collisions decrease the filling factor by a small amount. Thus, the blue curve closely follows the black line. In the low, intermediate, and highmass clouds, silica pebbles are significantly compressed to aggregates with volumefilling factors in the range of ⟨ φ ⟩ _{V} ≈ 0.22−0.43, regardless of their initial filling factor. Compression increases with increasing cloud mass, but the pebbles approach the maximum compression of φ_{max} ≈ 0.46 (Weidling et al. 2009; Güttler et al. 2010) only in the intermediate and highmass clouds. The maximum value is reduced by a factor 0.79 compared to φ_{2} in Table 1 because only the outer rim of the aggregate is compressed.
Ice pebbles, in contrast, show a different behaviour (see Fig. 3b). As in the case of silica, compression of ice pebbles in the verylow mass cloud is negligible. However, even here a few sticking collisions reduce the filling factor slightly. In the lowmass cloud, pebbles are compressed to ⟨ φ ⟩ _{V} ≈ 0.11. In simulation runs starting with φ_{0} ≳ 0.11, the ice pebbles retain their initial filling factor and the yellow line in Fig. 3b follows the black line. The same holds true for the intermediate and highmass clouds, but here the pebbles have at least ⟨ φ ⟩ _{V} ≈ 0.16, or approach their maximum filling factor of ⟨ φ ⟩ _{V} ≈ 0.23, according to the compression curve (see Fig. 1). Thus, φ = 0.11, φ = 0.16, and φ = 0.23 are lower boundaries for the compression of the ice pebbles.
Fig. 5 Packing fraction of pebbles as a function of the dusttoice ratio. The four panels correspond to different initial filling factors for the pebbles: a) of φ_{0} = 0.1; b) φ_{0} = 0.2; c) φ_{0} = 0.3; and d) φ_{0} = 0.4. Each panel shows the results for the four different cloud masses as solid curves: very lowmass (blue); lowmass (yellow); intermediatemass (green); and highmass (red). In panels a) and b), the curves for the very low and lowmass clouds are outside the plotting range, because the packing fraction becomes unphysically high (φ_{P} ≳ 1). In panel d), the green, blue, and yellow curves are overlapping because pebbles are barely compressed in these cases and thus have the same filling factor (see Fig. 4). The shaded area (violet) indicates the range of φ_{P} given by the random loose packing (RLP, φ_{P} = 0.55) and the random close packing (RCP, φ_{P} = 0.64). The error bars show the uncertainty of φ_{P}, based on the uncertainty of the volumefilling factor of the pebbles. Initial conditions, i.e. dusttoice ratio and cloud mass, which are consistent with a cometary bulk density of ρ_{c} = 0.5 g cm^{3}, result in a value for φ_{P} within the shaded area. Very low and low mass clouds require initially compact pebbles (φ_{0} = 0.4) and dusttoice ratios 3 ≲ ξ ≲ 7, whereas intermediate and highmass clouds need 3 ≲ ξ ≲ 9, regardless of initial filling factor. 
4.2.2. Mixed pebbles
Figure 4 shows the compression of mixed pebbles with different dusttoice ratios in the four clouds. In the very lowmass cloud, pebbles retain their initial filling factor regardless of dusttoice ratio. This is consistent with the findings for the pure ice and silica pebbles, which correspond to dusttoice ratios equal to 0 and ∞, respectively, in the very lowmass cloud (see Fig. 3). Pebble compression increases as the cloud mass increases from low to high. Additionally, increasing the dusttoice ratio also leads to an increased compression, because a higher dust content makes the pebbles weaker. The final filling factors are generally higher than for pure ice, but lower than for pure dust. In the intermediate and highmass clouds, the final filling factors are in the range ⟨ φ ⟩ _{V} ≈ 0.21−0.40, but only pebbles with a dusttoice ratio of ξ ≳ 3−5 obtain a filling factor close to ⟨ φ ⟩ _{V} ≈ 0.40.
4.3. Dusttoice ratio
We relate the dusttoice ratio of the pebbles to the dusttoice ratio of the cometesimal and thereby constrain the initial mass of the pebble cloud in which the cometesimal has formed. This is possible because we assume that the cloud collapse results in a single object without pebbles being lost. Consequently, the dusttoice ratio of the final cometesimal is the same as the dusttoice ratio of a single pebble.
The bulk density of the cometesimal is (14)where is the compact density of the pebble, ⟨ φ ⟩ _{V} is the final filling factor of the pebble, and φ_{P} is the packing fraction of pebbles within in the cometesimal. According to Skorov & Blum (2012), the packing fraction is φ_{P} ≈ 0.6, which falls in the range between random loose packing (RLP) with φ_{P} = 0.55 and random close packing (RCP) with φ_{P} = 0.64.
In our simulations, the bulk density of the cometesimal is a parameter and set to ρ_{c} = 0.5 g cm^{3} (e.g. Sierks et al. 2015). We use Eq. (14)to calculate φ_{P} for the given compact density and final filling factor of the pebbles. The resulting value of φ_{P} is then compared to the range given by RLP and RCP. This constrains the initial cloud mass because the final filling factor of the pebbles depends on the pebble cloud mass. Thus, initial conditions, i.e. ξ and M, which are consistent with ρ_{c}, are supposed to have φ_{p} in the given range.
In Fig. 5, we plot the pebble packing fraction versus dusttoice ratio. In each panel, the solid lines connect the φ_{P}s calculated from the pebblefilling factors that were obtained in our simulations (see Table 2). The error bars take the uncertainty of the volumefilling factor of the pebbles into account. Different line colours are for the different cloud masses. The shaded area highlights the range of φ_{P} given by the limits of RLP and RCP. From Figs. 5a and b, it is evident that initially porous pebbles in a very low or lowmass cloud are not sufficiently compressed to be consistent with a cometary bulk density of ρ_{c} = 0.5 g cm^{3} since they would require an unphysically dense packing (≳1). Thus, the lines are out of range and not shown. An initial filling factor of 0.3, however, requires a dusttoice ratio that is significantly higher than 10 and dense packing for the very low and lowmass clouds (Fig. 5c). If the pebbles are initially compact (φ_{0} ≳ 0.4), then very low and lowmass clouds are also capable of reproducing the expected φ_{P} for 3 ≲ ξ ≲ 7 (see Fig. 5d). The intermediate and highmass clouds lead to sufficiently compressed pebbles, regardless of their initial filling factor (see Figs. 3 and 4 and Table 2). However, only pebbles with a dusttoice ratio in the range 3 ≲ ξ ≲ 9 fall into the shaded region, which is consistent with a cometary bulk density of 0.5 g cm^{3}. Although it is possible to constrain the dusttoice ratio, the cloud mass is still ambiguous: lowmass clouds and initially compact pebbles work as well as highmass clouds with any initial pebble porosity.
5. Discussion
In this section we discuss the results presented in the previous section and we outline the caveats of our model.
5.1. Simulation results
5.1.1. Pebble compression
Bouncing collisions account for almost 100% of the collisional outcomes, regardless of the initial conditions (cloud mass and pebble filling factor) and material (see Sect. 4.1). Collisional compression, therefore, determines the final value of ⟨ φ ⟩ _{V}. This can be explained by considering the typical collision velocity of pebbles, v, during collapse, which determines the compression pressure p ∝ ρ_{•}v^{2} (Güttler et al. 2010). The virial velocity, which is a proxy for the velocity dispersion of the pebbles – and hence the typical collision speeds – is, by definition, , where R is the instantaneous cloud radius, and M is the (constant) cloud mass. Initially, the cloud radius equals the Hill radius, R_{0} = R_{H} = a(M/ 3M_{⊙})^{1 / 3}. We consider the cometesimal to be formed in the EdgeworthKuiper belt at a distance of a = 40 au. The virial velocity at this heliocentric distance is then given by (15)where r is the radius of the cometesimal in km, and η = R/R_{0} is the normalised size of the cloud (see also Wahlberg Jansson & Johansen 2014).
From Eq. (15), we infer that the typical collision velocities are of the order of cm s^{1} because . Thus, we mostly expect bouncing collisions to occur (see Fig. 7 of Windmark et al. 2012a). If we also take the higher threshold velocities for sticking, bouncing, and fragmentation of ice (a factor 10 compared to silica, e.g. Gundlach & Blum 2015) into account, bouncing is the only relevant collision type because pebbles collide neither fast enough for fragmentation, nor slow enough for sticking, which agrees with our findings in Sect. 4.1.
As a consequence of frequent bouncing collisions, silica pebbles are compressed during the collapse towards a filling factor of φ ≈ 0.40 (Weidling et al. 2009; Güttler et al. 2010). Thus, the silica pebbles are compact at the end of the collapse, even if they are porous at the beginning, and the cometary material consists of silica pebbles with a mean filling factor of ⟨ φ ⟩ _{V} ≈ 0.40 (see Fig. 3), which is in agreement with Skorov & Blum (2012); Blum et al. (2014). Ice pebbles, however, are harder to compress than silica pebbles and their maximum filling factor is lower (see Fig. 1). Therefore, for the same cloud mass, ice pebbles retain a higher porosity than silica pebbles, in agreement with Fig. 3. Increasing the cloud mass enhances the number of compact ice pebbles, because the velocity dispersion of pebbles, and thus the compression pressure, increases, as can be seen from Eq. (15).
Ice pebbles keep their low filling factor throughout the collapse, unless they had been compressed already before (in the very low, low and intermediate mass cloud cases), or the cloud mass, M, is sufficiently high (see Fig. 3b). However, 500 km is approximately the size of Ceres and too large for a comet, unless the cloud fragments during the collapse, which produces a size distribution of objects. Comets would then represent the lowmass tail of this distribution. We intend to address this scenario in future work.
Although the ice pebbles remain more porous than the silica pebbles, they do not retain an arbitrarily low filling factor. The lower threshold of φ is ⟨ φ ⟩ _{V} ≈ 0.11 for the lowmass cloud and ⟨ φ ⟩ _{V} ≈ 0.16 for the intermediatemass cloud. This is plausible because the pressure determines the compression in a bouncing collision. For fixed pressure, a pebble can only be compressed if its filling factor is below the compression curve (see Fig. 1). Otherwise, the pebble is too compact and further compression requires a higher pressure. As argued above, higher compression pressure requires a higher cloud mass. Thus, only pebbles with a filling factor less than a certain threshold can be compressed, which explains why the lower boundaries given above shift to higher values for increasing cloud mass, as seen in Fig. 3b. Mixed pebbles fall in between the limiting cases of pure ice and pure silica (see Fig. 4).
5.1.2. Dusttoice ratio and comet density
Based on the assumptions of mixed pebbles and a cloud collapse into a single cometesimal, a dusttoice ratio in the range of between 3 ≲ ξ ≲ 9 is necessary to be consistent with a cometary bulk density of ρ_{c} = 0.5 g cm^{3} and a pebble packing of φ_{P} ≈ 0.6.
From observations of cometary dust trails by the Infrared Astronomical Satellite (IRAS), Sykes & Walker (1992) inferred a dusttovolatile ratio of ~3. The comet 9P/Tempel 1 was visited by the NASA Deep Impact spacecraft and an impactor, launched from the spacecraft, hit the surface on 4 July 2005. From the ejected material, consisting of dust and water ice, Küppers et al. (2005) estimated a dusttoice ratio of >1 for the nucleus of 9P/Tempel 1. However, the ejecta originated from a nearsurface and ice depleted layer (the result of outgassing most likely), which was not representative for the bulk material of the comet. Thus, the ice proportion within the nucleus might still be higher, lowering the total dusttoice ratio. Based on GIADA (Grain Impact Analyser and Dust Accumulator) measurements aboard the Rosetta spacecraft, Rotundi et al. (2015) have recently estimated the dusttogas ratio of the comet 67P/ChuryumovGerasimenko as being 6 ± 2 for water ice only, and 4 ± 2 when CO and CO_{2} are taken into account. Our simulations are comparable with those measurements if comets form in pebble clouds of intermediate or highmass or in very low and lowmass clouds, but with initially compact pebbles. This is because, in those cases, we find a dusttoice ratio of >1, which falls into the range 1 <ξ< 10 for a typical cometary bulk density (see Fig. 5). However, we can exclude combinations of very low and lowmass clouds with fluffy pebbles as the initial conditions for cometesimal formation, because dusttoice ratios of ≫10 are not observed in comets. Bodies with such high dusttoice ratios would resemble rocks or maybe rocky asteroids with only a very little amount of ice.
With our simulations of either only ice or only silica pebbles we also cover two extreme cases of an ice and a silica cometesimal, respectively. The cloud contracts until it reaches a density of 0.5 g cm^{3}. Thus, we can then check whether the values obtained for the mean filling factor of the pebbles allow for a cometesimal entirely composed of only one of the two materials by calculating the packing fraction according to (14)for ξ = 0 and ξ = ∞, respectively. In our simulations, the first case of only ice always produces φ_{P} ≳ 1. This excludes pure ice cometesimals because the typical comet density of 0.5 g cm^{3} cannot be achieved with a packing fraction of φ_{p} ≈ 0.6 for the pebbles inside the comet (Skorov & Blum 2012; Blum et al. 2014, see Fig. 3a). The second case with only silica pebbles yields packing fractions that are either too high (≳0.7), or too low (≲0.5). However, in the very low and lowmass clouds, φ_{P} ≈ 0.6 is possible for 0.15 ≲ φ_{0} ≲ 0.3. From a physical perspective, a comet consisting of only dust is unlikely, because volatiles like water ice are essential, e.g. to explain activity (see Blum et al. 2014, 2015). On the other hand, only ice is also unlikely, because observations clearly show a certain amount of dusty material present in comets, e.g. dust trails or grains in the coma (Sykes & Walker 1992; Rotundi et al. 2015). Furthermore, it can be argued that a comet density that is different to our choice of 0.5 g cm^{3} is also covered by recent density estimates of comets (e.g. Blum et al. 2006; A’Hearn 2011; Sierks et al. 2015) because all methods rely on certain assumptions, e.g. density estimates based on orbital changes that are due to nongravitational forces require a detailed model of outgassing of volatiles. However, the density measurement of ρ_{c} = 0.470 ± 0.045 g cm^{3} (Sierks et al. 2015) for the Jupiter family comet 67P/ChuryumovGerasimenko by the Rosetta spacecraft points towards a density near the value that is assumed here. Thus, our results imply that a mixture of ice and dust is necessary to match the observations.
5.2. Cloud mass, pebble filling factor, and comet size
Although our model constrains the dusttoice ratio required to match observations to a range between 3 and 9, the initial cloud mass and the initial pebble filling factor are ambiguous. Intermediate and highmass clouds, regardless of the initial pebble filling factor work as well as very low and lowmass clouds with initially compact pebbles (see Sect. 4.3). Therefore, it is necessary to discuss this result with respect to observed comet sizes and pebble formation in more detail. We know from observations that cometary nuclei typically have sizes ranging from 1 km to 10 km (Lamy et al. 2004; A’Hearn 2011). This size range is consistent with our model of cometesimal formation in very low and lowmass clouds, which correspond to objects of size 0.5 km and 5 km, assuming a bulk density of 0.5 g cm^{3}, respectively. Thus, the cloud collapse directly leads to a comet. However, this requires initially compact pebbles. Starting with μmsized silica monomers, the growing aggregates encounter a bouncing barrier, which leads to the required compression (Zsom et al. 2010). However, based on the fact that ice is stickier, harder to fragment, and harder to compress (Gundlach et al. 2011; Gundlach & Blum 2015; Hill et al. 2015), Okuzumi et al. (2012) showed that 0.1 μmsized ice monomers grow to very porous aggregates, which do not experience significant collisional compression. Thus, we leave the issue of whether or not compact icy aggregates can form for future studies.
Our simulations also show that intermediate and highmass clouds are capable of reproducing the observed properties of cometary nuclei. The corresponding sizes of the cometesimals are 50 km and 500 km, respectively. In these cases, the pebbles are not required to be initially compact as they are compressed during the collapse. However, compared to the observed comet sizes of 1 km to 10 km, these cometesimals are too big. Thus, the cometesimal needs to fragment into several nuclei. One possibility enabling this to happen is during the collapse of the cloud by fragmentation into subclumps. The smaller subclumps then collapse into comet nuclei. Since the resulting objects are bound to the same cloud, binary formation should be possible (Nesvorný et al. 2010). An alternative possibility is the fragmentation of cometesimals in mutual collisions of big cometesimals after the collapse (Morbidelli & Rickman 2015). However, cometesimal collisions which are energetic enough to disrupt the colliding bodies can also be expected to significantly change the material properties, such as the density. As a consequence, the small fragments would neither resemble comets, nor would they carry a lot of information about the formation process of the cometesimal.
In summary, for this work on comet formation: in the framework of the streaming instability followed by the gravitational collapse of pebble clouds, a dusttoice ratio exceeding unity is required to match the typical cometary bulk density of 0.5 g cm^{3}, owing to the compression behaviour of the ice/dustmixed pebbles during the collapse. This indicates the formation of cometesimals in intermediate to highmass clouds. However, this implies that either the cloud or the cometesimal needs to fragment to form kmsized comet nuclei. Lower mass clouds are possible, if the pebbles already start off compact. However, very low and lowmass clouds in combination with fluffy pebbles can be excluded as initial conditions for cometesimal formation because, in these cases, the pebble packing becomes unphysically high (φ_{P} ≳ 1). Considering the conditions in the protoplantetary disk, which probably produces compact mixed pebbles with a filling factor in the range φ ~ 0.4 for pure silica and φ ~ 0.27 for pure ice (see Fig. 1), our model predicts a dusttoice ratio of ξ ~ 3−9 for comets. However, this needs to be investigated in more detail in future studies because it is not clear whether or not ice aggregates encounter a bouncing barrier in the way that silica particles do, while growing in the protoplanetary disk.
5.3. Caveats
We are aware of a number of caveats in our approach and try to draw some attention to these issues.
5.3.1. Pebble velocity distribution
First of all, as is done in Wahlberg Jansson & Johansen (2014), we assume that the pebble speeds follow a Maxwellian velocity distribution. This is only true if the collisions are frequent and fully elastic (ideal gas). However, the key ingredient of the collapse model is the inelasticity of pebble collisions. Sticking collisions are fully inelastic, whereas bouncing collisions can be either fully elastic or dissipate a certain fraction of the collision energy, which is characterised by the coefficient of restitution. Therefore, the velocity distribution function cannot be Maxwellian. The pebbles behave like a granular gas, which is an ensemble of macroscopic particles evolving under inelastic collisions. It can be shown analytically that, for a homogeneous granular gas without external forces, the velocity distribution function is close to Maxwellian (van Noije & Ernst 1998). The peak shifts slightly to lower velocities and the highvelocity tail is overpopulated because inelastic collisions tend to cool the system (van Noije & Ernst 1998). However, gravity is important in the collapsing cloud and the density distribution of pebbles is not necessarily homogeneous, e.g. as a result of clustering (Goldhirsch & Zanetti 1993). Bearing this in mind, more work on the velocity distribution of pebbles has to be carried out before the assumption of a Maxwellian velocity distribution is revisited.
5.3.2. Pebble size and triggering the streaming instability
Another uncertainty of the model is the pebble size we use. The streaming instability is sensitive to the local metallicity Z (i.e. the ratio of the dust and the gas column densities) of the protoplanetary disk and to the Stokes number of the pebbles. In simulations that do not take stratification of the disk that is due to the vertical component of the Sun’s gravity into account, particles with St ~ 1 are concentrated most effectively. On the other hand, in stratified simulations, pebbles with St ~ 0.3 show a strong clustering behaviour (Johansen et al. 2007, 2014). However, Bai & Stone (2010) reported clustering of pebbles via the streaming instability for St ≳ 0.01, if Z was significantly higher than the nominal value of the minimum mass solar nebula, Z_{MMSN} = 0.01 (Hayashi 1981).
The Stokes number depends on the properties of the gas and the pebbles. A minimum mass solar nebula model applied to 40 au shows that pebbles with radius s ≲ 800 m reside in the Epstein drag regime. The Stokes number can thus be written as (16)where ρ_{•} is the material density of the pebble, s is the size in cm, and a is the semimajor axis in au. Using the volumefilling factor, φ, and indicating compact values with a ^{∗}, we obtain (17)because ρ_{•} ∝ φ, and s ∝ φ^{− 1 / 3} by definition of the filling factor. In our simulations, we keep , s^{∗}, and a fixed, but vary φ, which means St ∝ 0.7 × φ^{2 / 3}. We therefore cover a range in Stokes numbers with our simulations and for φ ≲ 0.001, both ice and silica pebbles have St ≲ 0.01, which is lower than the critical value found in Bai & Stone (2010). Those pebbles would thus not trigger the streaming instability. A Stokes number of 0.3 (Johansen et al. 2007, 2014), on the other hand, corresponds to a filling factor of ~0.28 for silica. Thus, cmsized silica pebbles can trigger the SI, if they are compact as predicted by Zsom et al. (2010). To include the possibility of very fluffy ice aggregates, as mentioned in Okuzumi et al. (2012) and Kataoka et al. (2013), we have to take into account that those aggregates have to be very large to obtain a suitable Stokes number for streaming instability, e.g. St = 0.3 and φ = 0.001 requires msized ice pebbles. Compact ice pebbles with φ = 0.4 of ~cm in size, on the other hand, can also trigger the streaming instability. The same issues also apply for mixed pebbles, because their densities fall in the range between pure ice (1 g cm^{3}) and pure silica (3 g cm^{3}). This problem of initial conditions has to be managed more rigorously in future studies.
5.3.3. Selfgravity of the cometesimal
In our model, pebbles are collisionally compressed during the collapse. Once the comet is formed, compression by selfgravity becomes significant. Assuming a comet of constant density in hydrostatic equilibrium, the central pressure can be estimated as (18)where G is the gravitational constant, ρ_{c} is the comet density, and r is the comet radius (Blum et al. 2006). For a typical comet density of ρ_{c} = 0.5 g cm^{3}, a 0.5 kmsized object has a central pressure of P_{c} ≈ 9 × 10^{1} dyn cm^{2}, a 5 kmsized objects has P_{c} ≈ 9 × 10^{3} dyn cm^{2}, a 50 kmsized objects already has P_{c} ≈ 9 × 10^{5} dyn cm^{2}, and a 500 kmsized body has an even higher central pressure of 9 × 10^{7} dyn cm^{2}. Comparing these values to the compression curves for ice and silica pebbles shown in Fig. 1, we see that gravitational compression indeed plays a role for objects r ≳ 5 km. Consequently, the filling factors found for the pebbles must be regarded as lower limits. This also includes the mixed pebbles, because the pure silica and ice cases are limiting cases for dusttoice ratios of 0 and ∞, respectively.
Furthermore, we implicitly assume that the cometesimals are primordial objects. Further evolution of these icy bodies, like collisions between cometesimals or mass loss owing to activity, is excluded in our analysis. Although, those processes would arguably only affect the upper layer of the body.
5.3.4. Collision model for mixed pebbles
We constructed a collision model for mixed pebbles by interpolating physical properties, such as sticking, fragmentation, or compression of the pebbles between the two cases of pure ice and pure silica. The interpolation scheme is based on the fractional abundance of silica monomers within the pebble, because the collision outcome and the compression depends on the ability of the collision to rearrange or break monomer bonds (see Sect. 3.3.2 and Dominik & Tielens 1997; Blum & Wurm 2000). Thus, the approach is suited for pebbles that consist of a homogeneous mixture of μmsized silica and ice monomers. This is a reasonable assumption, because of the stochastic coagulation process in the solar nebula. In contrast, pebbles with a refractory core and an icy mantle are assumed to behave like ice pebbles but with enhanced mass, which is mainly determined by the silica core. However, we do not consider coremantle pebbles here.
As described in Sect. 3.3.2, our approach does not take compositional changes of the pebble into account. This simplification breaks down if the collisions become energetic enough to grind down pebbles into monomers. The following accretion of single monomers, which are either silica or ice, leads to aggregates with a different dusttoice ratio to the progenitor pebbles. However, the collision velocities in the collapsing pebble cloud are not high enough for this to be significant. In agreement with Wahlberg Jansson & Johansen (2014), we observe collision velocities, which are close to the virial velocity, in the less massive clouds. The velocities are of the order of a few cm s^{1} (see Eq. (15)), which lead to bouncing collisions only. In more massive clouds, the velocities are initially higher due to the higher virial velocity, which are 1−2 m s^{1}, see Eq. (15). However, the rapid energy dissipation, which is due to the onset of fragmentation in this velocity range (see Eqs. (5)−(7)), leads to velocities significantly lower than virial velocity again and also lower than the fragmentation threshold. Thus, pebbles and small fragments again undergo bouncing collisions. For decreasing dusttoice ratios, i.e. increasing ice content, the threshold velocities shift to even higher values and catastrophic fragmentation of pebbles is even less likely. However, for other applications, where velocities are higher and a significant amount of monomers is produced, a recipe for changing the composition of the aggregate might be necessary.
Furthermore, the compression of the pebbles also has to be performed with care. We use the omnidirectional dynamical compression curve for silica pebbles given by Güttler et al. (2010; see dotted line in Fig. 1). However, for the ice aggregates, it remains inconclusive whether the compression measured in the laboratory experiments is unidirectional or omnidirectional (see Sect. 2). Furthermore, the prescription for the porosity evolution given in Weidling et al. (2009) and Güttler et al. (2010) is tailored to fit silica particles. Adapting this model for ice by exchanging the compression curve only, is clearly not satisfactory. The same applies to mixed pebbles, where we combined the collision model of silica and ice. Therefore, we need more laboratory experiments for developing a model of ice compression in bouncing collisions, which is as detailed as the silica model.
5.3.5. Coefficient of restitution of pebbles
In our simulations, we chose a value of ~0.7 for the coefficient of restitution, which is suitable for dust aggregates and grazing collisions (Blum & Münch 1993; Weidling et al. 2012). Weidling & Blum (2015) performed microgravity experiments with mmsized dust aggregates and found that a coefficient of restitution uniformly distributed in a range between 0.29 and 0.81 with a mean value of 0.55 explains their experiments best. On the other hand, microgravity experiments with mm to cmsized ice pebbles have shown that the values for the coefficient of restitution are uniformly distributed between 0.06 and 0.84 with a mean value of 0.45 (Heißelmann et al. 2010) or between 0.08 and 0.65 with a mean value of 0.36 (Hill et al. 2015). Thus, our choice of coefficient of restitution is in agreement with both types of materials.
However, changing the coefficient of restitution slightly affects the final filling factor of the pebbles. A lower value increases the amount of energy dissipated in a single collision. As a consequence, the cloud collapse proceeds faster, the collision velocities are lower, and also the total number of collisions is lower. Thus, the pebbles are compressed less. From Eq. (14)we see that a lower filling factor of the pebbles increases the necessary packing fraction of the pebbles in the comet to obtain a comet density of 0.5 g cm^{3}. Hence, the curves in Fig. 5 shift slightly upwards and the range of dusttoice ratios shifts to higher values. However, this effect has no significant influence on our results for the very low and lowmass clouds, because the pebbles must be compact already at the beginning. In the intermediate and highmass clouds, the coefficient of restitution effect is partially compensated for by the higher collision velocities owing to the higher cloud mass, which leads to stronger compression. However, since pebble compression is determined by rearranging monomers, a lower coefficient of restitution should lead to a higher pebble compression. Our collision model does not capture this effect, because the pebble compression is determined independently from the coefficient of restitution, which only enters the collision velocity. Hence, we plan to investigate the influence of the coefficient of restitution in more detail in future work.
6. Conclusion
We investigated the collapse of a pebble cloud based on the model developed by Wahlberg Jansson & Johansen (2014) and followed the porosity evolution of the pebbles by using the collision model of Windmark et al. (2012a) for dust, which we adapted for ice by scaling the threshold velocities up by a factor 10 (see e.g. Gundlach & Blum 2015). Furthermore, a recipe for the evolution of the filling factor of dust pebbles in sticking and bouncing collisions from Güttler et al. (2010) was used. In addition, we developed an interpolated collision model for mixed pebbles containing both ice and dust. We also present novel results from laboratory experiments measuring the compressive strength of aggregates composed of μmsized water ice particles. These measurements are used to derive the evolution of the filling factor of the silica, ice, and mixed pebbles during the gravitational collapse.
We tested four initial cloud masses ranging from very low (2.6 × 10^{14} g) to high (2.6 × 10^{23} g), corresponding to cometesimals with radii varying from 0.5 km to 500 km and a bulk density of 0.5 g cm^{3}. We start with 1 cmsized pebbles that possess initial filling factors of between φ_{0} = 0.001 (very porous) and φ_{0} = 0.4 (compact), and investigated how the filling factor changes during the collapse, depending on the cloud mass. We conducted three sets of simulations, one with only silica pebbles, a second one with only ice pebbles, and a third one with mixed pebbles with varying dusttoice ratios. Additionally, we calculated the pebble packing required to produce the typical cometary bulk density of 0.5 g cm^{3} and investigate how this depends on the cloud mass, initial filling factor, and dusttoice ratio (see Fig. 5).
Our main findings can be summarised as follows:

The laboratory measurements of the compressive strength of granular ice show that water ice aggregates are harder to compress than silica aggregates (see Fig. 1). This results in a lower volumefilling factor for the ice aggregates compared with silica aggregates, if compressed with the same pressure.

Except for the very lowmass cloud, silica pebbles are always compressed and obtain volumefilling factors in the range of ⟨ φ ⟩ _{V} ≈ 0.22 − 0.43 at the end of the collapse, regardless of φ_{0} (see Fig. 3a).

Ice pebbles experience no significant compression in the very lowmass cloud. They are compressed to ⟨ φ ⟩ _{V} ≈ 0.11 and ⟨ φ ⟩ _{V} ≈ 0.17 in low and intermediatemass clouds, respectively; in highmass clouds, ice pebbles end up with ⟨ φ ⟩ _{V} ≈ 0.23 (see Fig. 3b). These values are lower boundaries for the pebble volumefilling factor. If the pebbles are more compact at the beginning, they retain their initial filling factor.

Mixed pebbles obtain filling factors in between the values for pure ice and pure silica.

Pebble compression, whereby the pebble packing is in agreement with the limits given by random loose packing (φ_{P} = 0.55) and random close packing (φ_{P} = 0.64) for a cometesimal with a bulk density of ρ_{c} = 0.5 g cm^{3}, can be achieved in two indistinguishable ways:

in very low and lowmass clouds for initially compact pebbles,

and in intermediate and highmass clouds, regardless of the initial filling factor of the pebbles.

The latter implies that one cometesimal forms more than one comet nucleus. In any case, a dusttoice ratio in the range 3 ≲ ξ ≲ 9 is necessary to match the observed bulk properties of comet nuclei.
Acknowledgments
We thank Karl Wahlberg Jansson for fruitful discussions. The ice compression experiments were performed under DLR grant 50WM1236. We thank Tobias Eckhardt for his help in the laboratory. We also thank the editor, the anonymous referees, and the language editor for their comments, which helped to improve the quality of this manuscript. S. Lorek is a member of the International MaxPlanck Research School for Solar System Science at the MaxPlanck Institute for Solar System Research, Göttingen.
References
 A’Hearn, M. F. 2011, ARA&A, 49, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Münch, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., Schräpler, R., Davidsson, B. J. R., & TrigoRodríguez, J. M. 2006, ApJ, 652, 1768 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2014, Icarus, 235, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2015, Icarus, 248, 135 [NASA ADS] [CrossRef] [Google Scholar]
 BockeléeMorvan, D., Gautier, D., Hersant, F., Huré, J.M., & Robert, F. 2002, A&A, 384, 1107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Connelly, J. N., Bizzarro, M., Krot, A. N., et al. 2012, Science, 338, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Cuzzi, J. N., & Weidenschilling, S. J. 2006, in ParticleGas Dynamics and Primary Accretion, eds. D. S. Lauretta, & H. Y. McSween (Tucson: University of Arizona Press), 353 [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Drążkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drążkowska, J.,Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldhirsch, I., & Zanetti, G. 1993, Phys. Rev. Lett., 70, 1619 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717 [NASA ADS] [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]
 Hanner, M. S., & Bradley, J. P. 2004, Composition and mineralogy of cometary dust, eds. M. C. Festou, H. U. Keller, & H. A. Weaver, 555 [Google Scholar]
 Hayashi, C. 1981, Progress of Theoretical Physics Suppl., 70, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Heißelmann, D., Blum, J., Fraser, H. J., & Wolling, K. 2010, Icarus, 206, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, C. R., Heißelmann, D., Blum, J., & Fraser, H. J. 2015, A&A, 573, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547 [Google Scholar]
 Jost, B., Gundlach, B., Pommerol, A., et al. 2013, Icarus, 225, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Küppers, M., Bertini, I., Fornasier, S., et al. 2005, Nature, 437, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Lamy, P. L., Toth, I., Fernandez, Y. R., & Weaver, H. A. 2004, The sizes, shapes, albedos, and colors of cometary nuclei, ed. G. W. Kronk, 223 [Google Scholar]
 McKeegan, K., Aléon, J., Bradley, J., et al. 2006, Science, 314, 1724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Morbidelli, A., & Rickman, H. 2015, A&A, 583, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rotundi, A., Sierks, H., Della Corte, V., et al. 2015, Science, 347, 3905 [Google Scholar]
 Schräpler, R., Blum, J., von Borstel, I., & Güttler, C. 2015, Icarus, 257, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, 1044 [NASA ADS] [CrossRef] [Google Scholar]
 Skorov, Y., & Blum, J. 2012, Icarus, 221, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Sunshine, J. M., Groussin, O., Schultz, P. H., et al. 2007, Icarus, 191, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Sykes, M. V., & Walker, R. G. 1992, Icarus, 95, 180 [NASA ADS] [CrossRef] [Google Scholar]
 van Noije, T., & Ernst, M. 1998, Granular Matter, 1, 57 [Google Scholar]
 Wahlberg Jansson, K., & Johansen, A. 2014, A&A, 570, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidling, R., & Blum, J. 2015, Icarus, 253, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Weidling, R., Güttler, C., Blum, J., & Brauer, F. 2009, ApJ, 696, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Weidling, R., Güttler, C., & Blum, J. 2012, Icarus, 218, 688 [NASA ADS] [CrossRef] [Google Scholar]
 Weissman, P. R., Asphaug, E., & Lowry, S. C. 2004, Structure and density of cometary nuclei, ed. G. W. Kronk, 337 [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012a, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012b, A&A, 544, L16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012c, A&A, 548, C1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [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]
All Tables
Volumeweighted mean filling factor of pebbles for all runs performed in this study.
All Figures
Fig. 1 Compression curves of the granular water ice samples consisting of micrometresized water ice particles (coloured curves). The initial volumefilling factor of the samples varies between 0.08 and 0.12. The dashed curve shows the best fit to the ice measurements. For comparison, the dotted line is the omnidirectional compression curve for silica particles used in the simulations (Güttler et al. 2009, 2010). 

In the text 
Fig. 2 Resolution study. Volumeweighted mean filling factor of pebbles as a function of the number of representative particles for different initial filling factors (φ_{0} = 0.01, 0.20, and 0.40) for the lowmass cloud. Increasing the number of particles does not change the final filling factor. 

In the text 
Fig. 3 Pebble compression for pure ice and silica pebbles. Volumeweighted mean filling factor at the end of collapse as a function of initial filling factor: a) silica pebbles; b) ice pebbles. The four cases of a very lowmass (blue), a lowmass (yellow), an intermediatemass (green), and highmass (red) cometesimal are shown. Along the black dotted line, the final filling factor equals the initial filling factor. The error bar indicate the standard deviation of the pebble ensemble. Silica pebbles are significantly compressed for M> 2.6 × 10^{20} g. The maximum filling factor ice pebbles can obtain by compression is considerably lower than for silica, φ ~ 0.23 compared to φ ~ 0.43. 

In the text 
Fig. 4 Pebble compression for mixed pebbles. Volumeweighted mean filling factor at the end of collapse as a function of initial filling factor: a) very lowmass cloud; b) lowmass cloud; c) intermediatemass cloud; d) highmass cloud. Mixed pebbles with different dusttoice ratios are shown: 0.5 (blue, icedominated), 1.0 (yellow), 3.0 (green, equal amount of dust and ice), 5.0 (red), 7.0 (violet), and 10 (brown, dustdominated). Along the black solid line, the final filling factor equals the initial filling factor. The error bar indicate the standard deviation of the pebble ensemble. As for pure silica or ice pebbles, compression increases with increasing cloud mass. Furthermore, compression increases with increasing dusttoice ratio. Pebbles with low dusttoice ratios behave like ice pebbles, while pebbles with high dusttoice ratios behave more like silica pebbles. 

In the text 
Fig. 5 Packing fraction of pebbles as a function of the dusttoice ratio. The four panels correspond to different initial filling factors for the pebbles: a) of φ_{0} = 0.1; b) φ_{0} = 0.2; c) φ_{0} = 0.3; and d) φ_{0} = 0.4. Each panel shows the results for the four different cloud masses as solid curves: very lowmass (blue); lowmass (yellow); intermediatemass (green); and highmass (red). In panels a) and b), the curves for the very low and lowmass clouds are outside the plotting range, because the packing fraction becomes unphysically high (φ_{P} ≳ 1). In panel d), the green, blue, and yellow curves are overlapping because pebbles are barely compressed in these cases and thus have the same filling factor (see Fig. 4). The shaded area (violet) indicates the range of φ_{P} given by the random loose packing (RLP, φ_{P} = 0.55) and the random close packing (RCP, φ_{P} = 0.64). The error bars show the uncertainty of φ_{P}, based on the uncertainty of the volumefilling factor of the pebbles. Initial conditions, i.e. dusttoice ratio and cloud mass, which are consistent with a cometary bulk density of ρ_{c} = 0.5 g cm^{3}, result in a value for φ_{P} within the shaded area. Very low and low mass clouds require initially compact pebbles (φ_{0} = 0.4) and dusttoice ratios 3 ≲ ξ ≲ 7, whereas intermediate and highmass clouds need 3 ≲ ξ ≲ 9, regardless of initial filling factor. 

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.