© 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/Churyumov-Gerasimenko 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 rubble-pile-like assemblage of dust and ice, held together by self-gravity rather than material strength (Weissman et al. 2004).

Growth from μm-sized dust and ice particles to km-sized comets is difficult to explain because the process is affected by several growth barriers: the bouncing barrier for mm- to cm-sized particles (Zsom et al. 2010), the fragmentation barrier for cm- to m-sized bodies (Blum & Wurm 2008; Güttler et al. 2010), and the radial drift barrier, which peaks for m-sized 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 mass-transfer collisions (Windmark et al. 2012a; Dra¸żkowska et al. 2013). It was found that a 100 m-sized 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 low-velocity sticking collisions and thus grow past the barriers (Windmark et al. 2012b,c). However, the dating of calcium-aluminium-rich 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 km-sized bodies via sweep-up growth is not fast enough, especially at large semi-major 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 two-fluid instability arises in the protoplanetary disk because dust and gas rotate at different velocities as a result of the sub-Keplerian speed of the gas and the non-perfect 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 Ceres-sized planetesimals (Johansen et al. 2007; Bai & Stone 2010; Dra¸żkowska & Dullemond 2014). The SI is triggered once the local dust-to-gas 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 collision-type 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 equal-sized 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 volume-filling 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. Hit-and-stick collisions of μm-sized dust grains form fractal aggregates with low-filling 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 volume-filling 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/ice-mixed 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 dust-to-ice 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 $\begin{array}{ccc}\mathit{\phi }\left(\mathit{p}\right)\mathrm{=}{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}\frac{{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{1}}}{\exp \frac{{\log }_{\mathrm{10}}\mathit{p}\mathrm{-}{\log }_{\mathrm{10}}{\mathit{p}}_{\mathrm{m}}}{\mathrm{\Delta }}\mathrm{+}\mathrm{1}}\mathrm{·}& & \end{array}$(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, φ₁ and φ₂ are the minimum and the maximum volume-filling factors of the material, respectively. Depending on the nature of the compression experiments, i.e. unidirectional or omnidirectional, the maximum volume-filling 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⁵ 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⁴ dyn cm^-2 by comparing their measurements to the results of smoothed-particle 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 power-law function, $\begin{array}{ccc}\mathit{\phi }\left(\mathit{p}\right)\mathrm{=}\frac{{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}{\mathit{\phi }}_{\mathrm{2}}}{\mathrm{2}}\mathrm{·}{\left(\frac{\mathit{p}}{{\mathit{p}}_{\mathrm{m}}}\right)}^{\frac{{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{1}}}{{\mathit{\phi }}_{\mathrm{2}}\mathrm{+}{\mathit{\phi }}_{\mathrm{1}}}\mathrm{·}\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }\ln \mathrm{10}}}\mathit{,}& & \end{array}$(2)to account for volume-filling 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).

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 volume-filling 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 micrometre-sized water ice particles with a mean radius of 1.45 μm. After sedimentation of the micrometre-sized water ice particles in the sedimentation chamber, a porous water ice aggregate with an initial volume-filling 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 micrometre-sized water ice particles (coloured curves). The initial volume-filling 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 set-up 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¹ 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 non-compressed volumes of the sample.

A total of four measurements are performed with this experimental set-up (solid curves in Fig. 1). The initial volume-filling factor varies between 0.08 and 0.12. Relatively low compressions (<10³ dyn cm^-2) do not lead to an increase in the volume-filling factor of the water ice samples. However, higher compressions are sufficient to compress the granular water ice samples from the initial volume-filling factor to a volume-filling factor of ~0.27 (for a compression of 10⁶ dyn cm^-2). The experiments stop at 10⁶ dyn cm^-2, because the existing experimental set-up (we adopted the experimental set-up 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 micrometre-sized 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 volume-filling 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 set-up 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⁶ dyn cm^-2, but more experiments with a modified set-up 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 non-compressed 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¹⁷ 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¹⁷ 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): $\begin{array}{ccc}& & \mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{d}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{s}}}\right)}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{18}}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\\ & & \mathrm{\Delta }{\mathit{v}}_{\mathrm{bounce}}^{\mathrm{d}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{b}}}\right)}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{18}}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\\ & & \mathrm{\Delta }{\mathit{v}}_{\mathrm{frag}}^{\mathrm{d}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{p}\mathit{,}\mathrm{t}}}{{\mathit{m}}_{\mathrm{1.0}}}\right)}^{-0.16}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}\end{array}$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⁷ g describes the onset of fragmentation. By replacing m_1.0 with m_0.5 = 9.49 × 10¹¹ 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: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{cms}}^{\mathrm{p}}& \mathrm{=}& \frac{\mathrm{\Delta }\mathit{v}}{\mathrm{1}\mathrm{+}{\mathit{m}}_{\mathrm{p}}\mathit{/}{\mathit{m}}_{\mathrm{t}}}\mathit{,}\\ {\mathit{v}}_{\mathrm{cms}}^{\mathrm{t}}& \mathrm{=}& \frac{\mathrm{\Delta }\mathit{v}}{\mathrm{1}\mathrm{+}{\mathit{m}}_{\mathrm{t}}\mathit{/}{\mathit{m}}_{\mathrm{p}}}\mathrm{·}\end{array}$For equal-mass 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 volume-filling 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. $\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{i}}\mathrm{=}\mathrm{10}\mathrm{\times }\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{d}}$, 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/ice-mixed pebbles

Using pure ice and pure silica pebbles is a major simplification. There is evidence for μm-sized 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ée-Morvan et al. 2002; Hanner & Bradley 2004; Cuzzi & Weidenschilling 2006). The coagulation timescale to grow from μm-sized grains to cm-sized pebbles is expected to be short (≲10³ yr, Zsom et al. 2010) compared to the radial mixing timescale (≳10⁴ yr, Bockelée-Morvan 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 comet-forming region. Thereby, the monomers coagulate to form cm-sized dust/ice-mixed 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 core-mantle 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_dm_0,d and m_i = N_im_0,i, respectively, and the pebble composition is expressed in terms of the dust-to-ice ratio, ξ ≡ m_d/m_i. The dust-to-ice 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 dust-to-ice ratio ξ. The pebble mass is the sum of the masses of all dust and ice monomers: $\begin{array}{ccc}\mathit{m}& & \mathrm{=}{\mathit{m}}_{\mathrm{d}}\mathrm{+}{\mathit{m}}_{\mathrm{i}}\mathrm{=}{\mathit{\rho }}_{\mathrm{d}}{\mathit{V}}_{\mathrm{d}}\mathrm{+}{\mathit{\rho }}_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}}\mathit{,}\end{array}$(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 $\begin{array}{ccc}\mathit{m}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{•}}\mathrm{(}{\mathit{V}}_{\mathrm{d}}\mathrm{+}{\mathit{V}}_{\mathrm{i}}\mathrm{)}}{\mathit{\phi }}\mathit{,}& & \end{array}$(9)where φ is the pebble-filling 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 ρ_•: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{•}}\mathrm{=}\mathit{\phi }\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{d}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{i}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\xi }\mathrm{)}}{\mathit{\xi }{\mathit{\rho }}_{\mathrm{i}}\mathrm{+}{\mathit{\rho }}_{\mathrm{d}}}\mathrm{=}\mathit{\phi }{\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\mathit{,}& & \end{array}$(10)where ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}$ is the compact density without any void space between the monomers. Thus, the density of the pebble depends on the dust-to-ice ratio. For ξ = 0 (no dust), ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\mathrm{=}{\mathit{\rho }}_{\mathrm{i}}$ and for ξ = ∞ (no ice) we get ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\mathrm{=}{\mathit{\rho }}_{\mathrm{d}}$, 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 $\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{mixed}}\mathrm{=}\mathit{x}\hspace{0.17em}\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{d}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{x}\mathrm{)}\hspace{0.17em}\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{i}}\mathrm{=}\mathrm{(}\mathrm{10}\mathrm{-}\mathrm{9}\mathit{x}\mathrm{)}\hspace{0.17em}\mathrm{\Delta }{\mathit{v}}_{\mathrm{stick}}^{\mathrm{d}}\mathit{.}$(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 dust-to-ice ratio of ξ, the remnant and the fragments will have the same value for the dust-to-ice 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 volume-filling 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 $\begin{array}{ccc}{\mathit{\phi }}^{\mathrm{mixed}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\mathit{x}\hspace{0.17em}{\mathit{\phi }}^{\mathrm{d}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{x}\mathrm{)}\hspace{0.17em}{\mathit{\phi }}^{\mathrm{i}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{,}& & \end{array}$(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 volume-filling 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₀, which ensures the correction term vanishes for collisions between two massive aggregates, depends on the monomer mass m₀. For mixed pebbles, however, silica and ice monomers are present and thus we replace m₀ by the reduced monomer mass ${\mathit{m}}_{\mathrm{0}}^{-1}\mathrm{}\mathrm{\to }{\mathit{m}}_{\mathrm{0}\mathit{,}\mathrm{d}}^{-1}\mathrm{+}{\mathit{m}}_{\mathrm{0}\mathit{,}\mathrm{i}}^{-1}$ in our implementation of the Ormel et al. (2007) porosity model.

3.4. Initial conditions

With this set-up (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, φ₀, of the pebbles and the pebble cloud mass, M. In the first set, we use water ice pebbles with a compact density of ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\mathrm{=}{\mathit{\rho }}_{\mathrm{i}}\mathrm{=}\mathrm{1}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{cm}}^{-3}$. In the second set, we have silica pebbles with a compact density of ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\mathrm{=}{\mathit{\rho }}_{\mathrm{d}}\mathrm{=}\mathrm{3}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{cm}}^{-3}$. In the third set, we simulate mixed pebbles with varying dust-to-ice 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 high-mass cloud or case, respectively. The exact masses of the clouds are 2.6 × 10¹⁴ g (very-low mass), 2.6 × 10¹⁷ g (low mass), 2.6 × 10²⁰ g (intermediate mass), and 2.6 × 10²³ 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 hit-and-stick 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 hit-and-stick 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 dust-to-ice ratios between 0.5 (ice-dominated) and 10 (dust-dominated). 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. Volume-weighted 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 low-mass 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 volume-filling factor of the pebbles. Figure 2 shows the volume-weighted mean filling factor of ice and silica pebbles for an increasing number of representative particles in the low-mass 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 volume-filling factor after ~2000 collisions. The mean number of collisions per representative particle also exceeds 2000 in the intermediate- and high-mass cloud for n = 100, which leads, together with the higher collision velocities in these cases, to a volume-filling 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 low-mass 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. Volume-filling factor

Fig. 3 Pebble compression for pure ice and silica pebbles. Volume-weighted 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 low-mass (blue), a low-mass (yellow), an intermediate-mass (green), and high-mass (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 × 1020 g. The maximum filling factor ice pebbles can obtain by compression is considerably lower than for silica, φ ~ 0.23 compared to φ ~ 0.43.

Fig. 4 Pebble compression for mixed pebbles. Volume-weighted mean filling factor at the end of collapse as a function of initial filling factor: a) very low-mass cloud; b) low-mass cloud; c) intermediate-mass cloud; d) high-mass cloud. Mixed pebbles with different dust-to-ice ratios are shown: 0.5 (blue, ice-dominated), 1.0 (yellow), 3.0 (green, equal amount of dust and ice), 5.0 (red), 7.0 (violet), and 10 (brown, dust-dominated). 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 dust-to-ice ratio. Pebbles with low dust-to-ice ratios behave like ice pebbles, while pebbles with high dust-to-ice ratios behave more like silica pebbles.

We investigate the final volume-weighted mean filling factor of pebbles: $\begin{array}{ccc}\mathrm{⟨}\mathit{\phi }{\mathrm{⟩}}_{\mathit{V}}\mathrm{=}\frac{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{N}}_{\mathit{i}}{\mathit{V}}_{\mathit{i}}{\mathit{\phi }}_{\mathit{i}}}{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{N}}_{\mathit{i}}{\mathit{V}}_{\mathit{i}}}\mathit{,}& & \end{array}$(13)where N_i = M/ (nm_i) is the number of particles a representative particle represents, ${\mathit{V}}_{\mathit{i}}\mathrm{=}{\mathit{V}}_{\mathit{i}}^{\mathrm{\ast }}\mathit{/}{\mathit{\phi }}_{\mathit{i}}$ is the volume of a pebble (${\mathit{V}}_{\mathit{i}}^{\mathrm{\ast }}$ is the compact volume, i.e. the volume without void space), and φ_i is the volume-filling 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 volume-filling 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 low-mass 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 high-mass clouds, silica pebbles are significantly compressed to aggregates with volume-filling 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 high-mass clouds. The maximum value is reduced by a factor 0.79 compared to φ₂ 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 very-low mass cloud is negligible. However, even here a few sticking collisions reduce the filling factor slightly. In the low-mass cloud, pebbles are compressed to ⟨ φ ⟩ _V ≈ 0.11. In simulation runs starting with φ₀ ≳ 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 high-mass 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 dust-to-ice 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 low-mass (blue); low-mass (yellow); intermediate-mass (green); and high-mass (red). In panels a) and b), the curves for the very low- and low-mass 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 volume-filling factor of the pebbles. Initial conditions, i.e. dust-to-ice 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 dust-to-ice ratios 3 ≲ ξ ≲ 7, whereas intermediate- and high-mass clouds need 3 ≲ ξ ≲ 9, regardless of initial filling factor.

4.2.2. Mixed pebbles

Figure 4 shows the compression of mixed pebbles with different dust-to-ice ratios in the four clouds. In the very low-mass cloud, pebbles retain their initial filling factor regardless of dust-to-ice ratio. This is consistent with the findings for the pure ice and silica pebbles, which correspond to dust-to-ice ratios equal to 0 and ∞, respectively, in the very low-mass cloud (see Fig. 3). Pebble compression increases as the cloud mass increases from low to high. Additionally, increasing the dust-to-ice 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 high-mass clouds, the final filling factors are in the range ⟨ φ ⟩ _V ≈ 0.21−0.40, but only pebbles with a dust-to-ice ratio of ξ ≳ 3−5 obtain a filling factor close to ⟨ φ ⟩ _V ≈ 0.40.

4.3. Dust-to-ice ratio

We relate the dust-to-ice ratio of the pebbles to the dust-to-ice 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 dust-to-ice ratio of the final cometesimal is the same as the dust-to-ice ratio of a single pebble.

The bulk density of the cometesimal is $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{c}}\mathrm{=}{\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}\hspace{0.17em}\mathrm{⟨}\mathit{\phi }{\mathrm{⟩}}_{\mathit{V}}\hspace{0.17em}{\mathit{\phi }}_{\mathrm{P}}\mathit{,}& & \end{array}$(14)where ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}$ 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 dust-to-ice ratio. In each panel, the solid lines connect the φ_Ps calculated from the pebble-filling factors that were obtained in our simulations (see Table 2). The error bars take the uncertainty of the volume-filling 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 low-mass 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 dust-to-ice ratio that is significantly higher than 10 and dense packing for the very low- and low-mass clouds (Fig. 5c). If the pebbles are initially compact (φ₀ ≳ 0.4), then very low- and low-mass clouds are also capable of reproducing the expected φ_P for 3 ≲ ξ ≲ 7 (see Fig. 5d). The intermediate- and high-mass 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 dust-to-ice 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 dust-to-ice ratio, the cloud mass is still ambiguous: low-mass clouds and initially compact pebbles work as well as high-mass 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² (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, $\mathit{\sigma }\mathrm{\approx }\sqrt{\mathit{GM}\mathit{/}\mathit{R}}$, where R is the instantaneous cloud radius, and M is the (constant) cloud mass. Initially, the cloud radius equals the Hill radius, R₀ = R_H = a(M/ 3M_⊙)^{1 / 3}. We consider the cometesimal to be formed in the Edgeworth-Kuiper belt at a distance of a = 40 au. The virial velocity at this heliocentric distance is then given by $\begin{array}{ccc}\mathit{\sigma }\mathrm{=}\mathrm{2.9}\hspace{0.17em}{\left(\frac{\mathit{a}}{\mathrm{40}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{{\mathit{\rho }}_{\mathrm{c}}}{\mathrm{0.5}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{km}}\right)\hspace{0.17em}{\mathit{\eta }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}& & \end{array}$(15)where r is the radius of the cometesimal in km, and η = R/R₀ 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 $\mathit{v}\mathrm{~}\sqrt{\mathrm{2}}\mathit{\sigma }$. 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 low-mass 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 low-mass cloud and ⟨ φ ⟩ _V ≈ 0.16 for the intermediate-mass 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. Dust-to-ice ratio and comet density

Based on the assumptions of mixed pebbles and a cloud collapse into a single cometesimal, a dust-to-ice 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 dust-to-volatile 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 dust-to-ice ratio of >1 for the nucleus of 9P/Tempel 1. However, the ejecta originated from a near-surface 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 dust-to-ice ratio. Based on GIADA (Grain Impact Analyser and Dust Accumulator) measurements aboard the Rosetta spacecraft, Rotundi et al. (2015) have recently estimated the dust-to-gas ratio of the comet 67P/Churyumov-Gerasimenko as being 6 ± 2 for water ice only, and 4 ± 2 when CO and CO₂ are taken into account. Our simulations are comparable with those measurements if comets form in pebble clouds of intermediate- or high-mass or in very low- and low-mass clouds, but with initially compact pebbles. This is because, in those cases, we find a dust-to-ice 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 low-mass clouds with fluffy pebbles as the initial conditions for cometesimal formation, because dust-to-ice ratios of ≫10 are not observed in comets. Bodies with such high dust-to-ice 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 low-mass clouds, φ_P ≈ 0.6 is possible for 0.15 ≲ φ₀ ≲ 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 non-gravitational 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/Churyumov-Gerasimenko 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 dust-to-ice 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 high-mass clouds, regardless of the initial pebble filling factor work as well as very low- and low-mass 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 low-mass 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 μm-sized 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 μm-sized 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 high-mass 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 sub-clumps. The smaller sub-clumps 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 dust-to-ice 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/dust-mixed pebbles during the collapse. This indicates the formation of cometesimals in intermediate- to high-mass clouds. However, this implies that either the cloud or the cometesimal needs to fragment to form km-sized comet nuclei. Lower mass clouds are possible, if the pebbles already start off compact. However, very low- and low-mass 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 dust-to-ice 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 high-velocity 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 $\begin{array}{ccc}\mathrm{St}\mathrm{=}\mathrm{0.7}\mathrm{\times }\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{•}}}{\mathrm{3}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)\hspace{0.17em}\left(\frac{\mathit{s}}{\mathrm{1}\hspace{0.17em}\mathrm{cm}}\right)\hspace{0.17em}{\left(\frac{\mathit{a}}{\mathrm{40}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(16)where ρ_• is the material density of the pebble, s is the size in cm, and a is the semi-major axis in au. Using the volume-filling factor, φ, and indicating compact values with a ^∗, we obtain $\begin{array}{ccc}\mathrm{St}\mathrm{=}\mathrm{0.7}\mathrm{\times }{\mathit{\phi }}^{\mathrm{2}\mathit{/}\mathrm{3}}\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}}{\mathrm{3}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)\hspace{0.17em}\left(\frac{{\mathit{s}}^{\mathrm{\ast }}}{\mathrm{1}\hspace{0.17em}\mathrm{cm}}\right)\hspace{0.17em}{\left(\frac{\mathit{a}}{\mathrm{40}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}& & \end{array}$(17)because ρ_• ∝ φ, and s ∝ φ^{− 1 / 3} by definition of the filling factor. In our simulations, we keep ${\mathit{\rho }}_{\mathrm{•}}^{\mathrm{\ast }}$, 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, cm-sized 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 m-sized 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. Self-gravity of the cometesimal

In our model, pebbles are collisionally compressed during the collapse. Once the comet is formed, compression by self-gravity becomes significant. Assuming a comet of constant density in hydrostatic equilibrium, the central pressure can be estimated as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{c}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{,}& & \end{array}$(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 km-sized object has a central pressure of P_c ≈ 9 × 10¹ dyn cm^-2, a 5 km-sized objects has P_c ≈ 9 × 10³ dyn cm^-2, a 50 km-sized objects already has P_c ≈ 9 × 10⁵ dyn cm^-2, and a 500 km-sized body has an even higher central pressure of 9 × 10⁷ 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 dust-to-ice 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 μm-sized 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 core-mantle 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 dust-to-ice 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 dust-to-ice 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 mm-sized 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 cm-sized 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 dust-to-ice ratios shifts to higher values. However, this effect has no significant influence on our results for the very low- and low-mass clouds, because the pebbles must be compact already at the beginning. In the intermediate- and high-mass 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 μm-sized 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¹⁴ g) to high (2.6 × 10²³ 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 cm-sized pebbles that possess initial filling factors of between φ₀ = 0.001 (very porous) and φ₀ = 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 dust-to-ice 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 dust-to-ice 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 volume-filling factor for the ice aggregates compared with silica aggregates, if compressed with the same pressure.

• Except for the very low-mass cloud, silica pebbles are always compressed and obtain volume-filling factors in the range of ⟨ φ ⟩ _V ≈ 0.22 − 0.43 at the end of the collapse, regardless of φ₀ (see Fig. 3a).

• Ice pebbles experience no significant compression in the very low-mass cloud. They are compressed to ⟨ φ ⟩ _V ≈ 0.11 and ⟨ φ ⟩ _V ≈ 0.17 in low- and intermediate-mass clouds, respectively; in high-mass clouds, ice pebbles end up with ⟨ φ ⟩ _V ≈ 0.23 (see Fig. 3b). These values are lower boundaries for the pebble volume-filling 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 low-mass clouds for initially compact pebbles,

  • and in intermediate- and high-mass 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 dust-to-ice 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 Max-Planck Research School for Solar System Science at the Max-Planck Institute for Solar System Research, Göttingen.

References

All Tables

All Figures

Fig. 1 Compression curves of the granular water ice samples consisting of micrometre-sized water ice particles (coloured curves). The initial volume-filling 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. Volume-weighted 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 low-mass 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. Volume-weighted 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 low-mass (blue), a low-mass (yellow), an intermediate-mass (green), and high-mass (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 × 1020 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. Volume-weighted mean filling factor at the end of collapse as a function of initial filling factor: a) very low-mass cloud; b) low-mass cloud; c) intermediate-mass cloud; d) high-mass cloud. Mixed pebbles with different dust-to-ice ratios are shown: 0.5 (blue, ice-dominated), 1.0 (yellow), 3.0 (green, equal amount of dust and ice), 5.0 (red), 7.0 (violet), and 10 (brown, dust-dominated). 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 dust-to-ice ratio. Pebbles with low dust-to-ice ratios behave like ice pebbles, while pebbles with high dust-to-ice ratios behave more like silica pebbles.

In the text

Fig. 5 Packing fraction of pebbles as a function of the dust-to-ice 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 low-mass (blue); low-mass (yellow); intermediate-mass (green); and high-mass (red). In panels a) and b), the curves for the very low- and low-mass 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 volume-filling factor of the pebbles. Initial conditions, i.e. dust-to-ice 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 dust-to-ice ratios 3 ≲ ξ ≲ 7, whereas intermediate- and high-mass clouds need 3 ≲ ξ ≲ 9, regardless of initial filling factor.

In the text