© ESO, 2014

1. Introduction

In the process of forming planets the initially μm-sized dust and ice particles in the primordial protoplanetary disk grow to planets of sizes up to ~10⁴ km. This is believed to occur in distinct stages where different physical growth processes dominate (Safronov 1969). The growth from dust to planetesimals occurs at first via direct sticking, but is later aided by the increase in the local particle density by concentration in the turbulent gas (Johansen et al. 2014). Small particles can stick together via surface forces, but when they reach mm and cm sizes, growth stops as collisions lead to bouncing (Zsom et al. 2010) and fragmentation (Brauer et al. 2008; Birnstiel et al. 2009). The two leading theories for the growth mechanism beyond this point are 1) hierarchical coagulation and 2) gravitational instability.

Growth by coagulation can proceed via mass transfer as small impactors deposit part of their mass when hitting a large target particle (Wurm et al. 2005). The resulting growth is nevertheless slow, resulting in the formation of 100 m scale planetesimals in the asteroid belt in 1 Myr (Windmark et al. 2012); and growth rates would be even lower in the Kuiper belt.

Planetesimal formation via gravitational instability is also not without problems. Gravitational instability of a region in an unperturbed protoplanetary disk can be easily stopped by Kelvin-Helmholz turbulence in the mid-plane layer (Goldreich & Ward 1973; Weidenschilling 1980). However, the presence of gas in the disk can cause particles to clump together through streaming instability (Youdin & Goodman 2005) and lead to high particle overdensities. Simulations of this mechanism in protoplanetary disks result in the formation of gravitationally bound pebble clouds with masses comparable to ~100 km solid planetesimals or larger (e.g. Johansen et al. 2009, 2012).

An advantage with gravitational instability is that, if the internal angular momentum of a gravitationally unstable pebble cloud is large, it naturally forms a binary planetesimal with mass ratio of ~1 (Nesvorný et al. 2010). The components form from the same material, so they will also have the same chemical composition and colour, in agreement with what is observed for binary Kuiper belt objects (Benecchi et al. 2009). The Kuiper belt as a whole has a wide colour distribution, so this suggests that the components in binary Kuiper belt objects formed together; a similar co-formation explanation is not possible in models of binary formation through three-body encounters (Goldreich et al. 2002) in the hierarchical coagulation formation mechanism.

Nesvorný et al. (2010) modelled the particle cloud with an N-body method with perfect sticking. The goal of this paper is to use more realistic collision outcomes based on laboratory experiments (Güttler et al. 2010) to model the collapse. Inelastic collisions lead to energy dissipation and contraction of the cloud, which in turn releases gravitational energy into random pebble motion. As a pebble cloud collapses, the relative speed between the pebbles can, depending on the mass of the cloud, reach values large enough to result in fragmentation. Low-mass clouds experience no fragmentation and the planetesimal will be formed from primordial pebbles. The pebbles in high-mass clouds, on the other hand, are ground down to dust during the collapse. Depending on the degree of fragmentation, the porosity and density of the resulting planetesimal will vary. In the case of Kuiper belt objects, measured densities range from below 1 g cm^-3 for smaller bodies (radii ≲200 km) up to ~2.5 g cm^-3 for Eris (radius ~1000 km; Fig. 3 Brown 2013). The density of a Kuiper belt object is not only determined by the porosity, but also by the composition. A large ice fraction results in a lower density, suggesting that more massive objects are both less porous and have a higher rock fraction.

Observations of comet nuclei in the solar system indicate densities of ~0.5 g cm^-3 (Weissman et al. 2004; Davidsson et al. 2007) due to high porosity combined with large ice fraction. High porosities can be either mesoscale (rubble piles) or small scale (pebble piles). A pebble pile could originate from the collapse of a low-mass pebble cloud where the collision speeds never reach the fragmentation limit. Studies of the tidal break-up of the comet Shoemaker-Levy 9 showed that its internal strength was as low as expected for a rubble (or pebble) pile (Asphaug & Benz 1996). Observations of the comet 67P/Churyumov-Gerasimenko (Bauer et al. 2012) show that mm-sized pebbles are released from the surface of the comet. Skorov & Blum (2012) model this pebble release process with a two-layered comet nucleus consisting of an ice-free surface layer residing on top of a mixture of refractory pebbles and ice aggregates. These observations agree with the idea of comets being piles of mm- and cm-sized pebbles. The interior structure of comets will also be investigated by the Rosetta mission which will put a lander on 67P/Churyumov-Gerasimenko in November 2014.

An additional motivation for this paper is the NASA mission New Horizons, which is on its way to the outer regions of our solar system. It is scheduled to arrive at the Pluto-Charon system in July 2015 (Stern 2008). After a fly-by it is planned to encounter and investigate one or two additional Kuiper belt objects. The New Horizons probe will measure the mass, binarity, surface composition, and temperature in such fly-bys. As no currently known Kuiper belt objects are within the range of New Horizons after the fly-by of Pluto a survey is ongoing to detect possible candidates (Buie et al. 2012).

We have developed a model for simulating the collapse of gravitationally bound pebble clouds. With this model we investigate e.g. the collapse time and particle size distribution as a function of planetesimal size. The paper is organized as follows. We present our zero-dimensional model of the cloud, the model of collisional outcomes, and the numerical scheme in Sect. 2. The simulations of clouds with different masses and initial pebble sizes are presented in Sect. 3 and the results are discussed in Sect. 4. In Appendix A we present an analytic description of the collapse of a low-mass pebble cloud only experiencing bouncing collisions.

2. Model

We start with a gravitationally bound cloud of pebbles formed for example through the streaming instability. As the pebbles move around they dissipate energy in inelastic collisions with each other leading to a contraction of the cloud. Since gravitationally bound systems have a negative heat capacity, the pebbles will increase their speed when the cloud loses energy. This together with the increased density as the cloud contracts will lead to higher collision rates and faster energy dissipation. The result is a runaway collapse. With some approximations (e.g. only bouncing collisions, equal particle size, and immediate virialization) this collapse can be derived analytically (see Appendix A). Including more realistic collision physics requires numerical solution.

2.1. Parametrization of the cloud size and energy

Our model is zero-dimensional, so it requires that the cloud is completely uniform and that there is no net rotation. This is of course not completely physical but it is a first step towards a more complete model. We also assume that the cloud is in, or strives to get into, virial equilibrium. The cloud is initially set in virial equilibrium so that the initial total energy, E₀, is (1)Here T₀ is the initial total kinetic energy (the sum of the kinetic energy of each individual particle with mass m_i and speed v_i), M is the total mass of the pebble cloud, v_vir,0 is the virial speed calculated from the initial kinetic energy, U₀ is the initial potential energy, G is the gravitational constant, and R₀ is the initial radius of the cloud. During the collapse we can follow the cloud’s evolution through three parameters: η_eq, η, and η_K. The parameter η_eq is defined as the initial energy divided by the current total energy E(t), (2)The parameter η is defined as the current size, R(t), of the cloud divided by the initial size, (3)Finally η_K, which we use to generate relative speeds when the cloud is not in virial equilibrium, is defined as the initial kinetic energy divided by the current kinetic energy, T(t), (4)The cloud is bound and contracting, so η_eq and η have values between 0 and 1, while η_K can have any value since kinetic energy is both dissipated in collisions and released in the gravitational collapse. These three parameters do not need to have the same value since the cloud does not virialize immediately after a collision. We assume that the cloud virializes on a timescale of the free-fall time of the cloud. At a certain size of the cloud, η, we can get the maximum possible change in η, δη_max, by assuming that the cloud has been free-falling from R₀. This gives (5)where t_ff,0 is the initial free-fall time. Here we have used the free-fall speed at R starting from R₀, (6)The desired change in η after a dissipative collision is δη = η_eq − η. If | δη | < | δη_max | the new value of η is equal to the equilibrium value, otherwise the new value is the old plus δη_max, since there is not enough time for the system to virialize in the time step δt.

Fig. 1 Collisional outcomes as a function of projectile radius, a, collision speed, Δv, and relative particle size, f ≡ m2/m1. Green regions indicate coagulation, yellow regions bouncing, and red regions fragmentation. Outcome regions taken from Güttler et al. (2010).

When it comes to updating η_K we need to take into account that kinetic energy is both dissipated in the collisions and released when the cloud contracts, so using Eqs. (1), (3), and (4)we obtain (7)where δE is the energy lost in the collision (see Sect. 2.2), δU is the gravitational energy released, and δη is the previously calculated change in η (could be δη_max). The first term in Eq. (7)decreases the kinetic energy (energy dissipation in collisions) while the second term increases the kinetic energy (release of gravitational potential in the collapse). If the virialization is too slow (i.e. if the next collision occurs before the gravitational energy has been released), δη_K>δη and the particles will achieve subvirial velocities.

2.2. Collisional outcomes

For the cloud to collapse it needs to dissipate a lot of kinetic energy. In our model this is done by inelastic collisions between the particles. We use the tabulated laboratory measurements of Güttler et al. (2010) to get the outcome of a collision between a target particle (subscript 1) and a projectile (subscript 2). The outcome depends on the relative speed, Δv, and relative size, f ≡ m₂/m₁, of the particles. The collisions are split into two groups:

• 1.

  two similar-sized particles colliding or a larger particle, 2, colliding with a small particle, 1: f ≥ 0.1;

• 2.

  a small particle, 2, running into (colliding with) a larger particle, 1: f< 0.1.

We do a simple interpretation of Fig. 11 in Güttler et al. (2010), for collisions between compact particles:

Here C denotes coagulation, B bouncing, and F fragmentation. The threshold speed for sticking, v_stick, is taken from Güttler et al. (2010), as (8)where μ = m₁m₂/ (m₁ + m₂) is the reduced mass, a₀ is the monomer radius, and F_roll is the rolling force of the monomers¹. The coagulation outcome of collisions with f< 0.1 and 1 m s^-1≤ Δv< 25 m s^-1 is explained by mass transfer or penetration (Wurm et al. 2005). A map of collision outcomes in (Δv,a)-space is shown in Fig. 1.

Next we determine the amount of energy dissipated, δE, in a collision between target particle 1 and projectile 2. In the case of bouncing we use conservation of momentum and the coefficient of restitution, C_R, to get the change in energy, (9)For simplicity we assume that all kinetic energy in the direction of relative motion is dissipated, C_R = 0. Particle properties like mass and radius remain constant in a bouncing collision. In the case of coagulation we use the same equation but the difference from bouncing is that the mass of the target increases to and hence the radius also changes.

A constant C_R = 0 is not physical since the collisions might not be completely inelastic. Higa et al. (1996) investigate the coefficient of restitution for collisions between ice pebbles. They find that C_R depends on the collision speed and for low values (Δv ≲ 0.35 m s^-1) the collisions are almost elastic (C_R ~ 0.9). For higher velocities, on the other hand, C_R decreases with increasing collision speed and reaches values close to zero at a few m s^-1. Compared to our simulations these experiments investigate ice pebbles while we consider silicates. A non-zero coefficient of restitution would decrease the amount of energy dissipated in each bouncing collision which, in turn, would increase the timescales of the collapses. Another implication is that the clouds would get more time to virialize and the collision speeds becomes a bit higher. However, the value of the coefficient of restitution would need to be relatively high to give a significant effect to the collapse since it is the square of C_R that appears in Eq. (9).

Fig. 2 In our model pebbles are composed of μm-sized monomer particles. In collisions resulting in partial fragmentation monomers are removed one by one to simulate erosion. The binding energy of the fragment, , is larger (closer to zero) than the binding energy of the original particle, Ebind, since the fragment consists of fewer monomers. The number of monomers released, Nmon, is equal the number of contact surfaces the collision energy can break divided by three contact surfaces per monomer, Nmon = Ecoll/(3Eroll).

The effect of fragmentation is a bit different. In a collision between two particles with mass m₁ and m₂ we have, in the centre-of-mass frame, an available collision energy of (10)Next we assume that all particles are composed of μm-sized monomers of mass m₀. The energy with which any contact surface of two monomers is held together is (Dominik & Tielens 1997; Blum & Wurm 2000) (11)With Eqs. (10), (11)we get the number of bonds that can be broken by the collision as N_cs ~ E_coll/E_roll. For a compact aggregate of monomers, each monomer has, on average, three contact surfaces (Zsom & Dullemond 2008) so the energy required to completely fragment a particle is (12)If E_coll ≥ E_frag the whole particle is destroyed and μm-sized monomers are produced. The energy dissipated is δE = −E_frag, as not all collision energy is needed to fragment the target. If, however, E_coll<E_frag some part of the old particle survives with N′ monomers and ~3N′ contact surfaces. For simplicity we assume that in the collision the monomers are removed one by one (see sketch in Fig. 2). The mass of the fragment becomes (13)In this case the energy dissipated is δE = −E_coll, all the collision energy goes into partially fragmenting the target.

To make the simulations more physically realistic we add an impact parameter, b, to all collisions. The value of b is randomly picked between 0 and the sum of the radii (a₁ and a₂) of the two particles with higher probability for larger impact parameters (larger area), (14)The effect of a non-zero impact parameter is that the effective relative speed decreases. This means that the outcome of the collision can change and the energy dissipated will be lower.

2.3. Representative particle approach

To model the collapse of a planetesimal mass cloud of mm- and cm-sized pebbles, formed e.g. through the streaming instability (Johansen et al. 2009, 2012), it is not possible to follow the evolution of every single particle. A Pluto-mass planetesimal split into cm-sized pebbles results in N ~ 10²⁴ particles. Therefore we use the representative particle approach of Zsom & Dullemond (2008) where the particles are grouped into a smaller number (manageable by a computer) of superparticles, each representing a large number of particles with identical properties. We repeat the main steps of the Zsom & Dullemond (2008) algorithm here.

The representative particle approach is a Monte Carlo model that uses the collision rates of all possible collision pairs to find which particles collide and the times between collisions. Initially a large number, N, of physical particles are defined. From these we randomly select a smaller number, most simulations are done with n = 1000, representative particles which we follow the evolution of. Each representative particle, i, has its own properties like mass or velocity. One can think of a representative particle as one member of a swarm of identical physical particles: if the properties of the representative particle changes the properties of all the particles within the swarm changes. In the simulation the total mass of each of these swarms remains constant (M_i = M/n) so the number of physical particles in a swarm will change.

To follow the evolution of the cloud we need to know with what rate a representative particle collides with any other particle. Assume that n is large enough so the distributions of the properties of the representative particles is a good representation of the true distributions of all N physical particles. We then get the collision rate between representative particle i and a physical particle with the same properties as representative particle k to be (15)where n_k is the number density of physical particles represented by representative particle k, σ_ik is the cross-section between two such particles, and Δv_ik is the relative speed. In our simulations we generate the relative speeds from the kinetic energy of the cloud, η_K, and assuming a maxwellian distribution of the relative speeds, (16)where is the one-dimensional velocity dispersion in the cloud and we have included the fact that the relative speeds are, on average, a factor of larger than the individual speeds. For strongly dissipative collisions we do not expect energy equipartition to occur and hence the small particles do not necessarily have much higher velocities than the large particles. The actual velocity differences are not straightforward to obtain (Barrat & Trizac 2002) and we make the assumption that the particle velocities for the dissipative system is independent of particle mass. From Eq. (15)we calculate the rate with which representative particle i collides with any particle (it can collide with particles in its own swarm as well), (17)and the total rate of collisions for any of the representative particles, (18)From this total rate and the assumption that the collision process is a Poisson process we get the exponentially distributed time until next collision, (19)where rnd(seed) is a random number uniformly distributed between 0 and 1. Next we choose which particles collide with two more random numbers. First the representative particle i with the probability (20)and then the physical particle k, given representative particle i, with the probability (21)Knowing which two particles collide we calculate the outcome and effects of such a collision. For this we use the equations in Sect. 2.2 changing the subscripts 1 and 2 to i and k, i.e. the representative particle is the target and the physical particle is the projectile. We also need to take into consideration that all equations in Sect. 2.2 are only for one collision. Therefore, when calculating the energy dissipated we need to multiply the equations with the number of physical particles in the swarm, . Another difference arises when we pick the new mass of the representative particle after a fragmenting collision. If we have complete fragmentation the new mass is simply the monomer mass m₀. If, on the other hand, we have partial fragmentation, we pick the new mass from a mass-weighted bimodal distribution defined by the value (Eq. (13)). We pick a random number uniformly distributed between 0 and m_i, if it is ≤ the new mass is otherwise the new mass is m₀. Basically we randomly select a monomer in the original particle and see where it ends up.

After a collision, say representative particle j collided with physical particle l, the properties of the representative particle has changed and the rate matrix, Eq. (15), needs to be updated. Representative particle l has not been involved in the collision and since we are only looking at the evolution of the representative particles the only elements that needs to updated are the ones with representative particle i = j and physical particle k = j (one row and one column in r_i(k)). The symmetry is restored when representative particle l collides with physical particle j. Here we also take advantage of the η-parametrization of the cloud. After every collision the size of the cloud and the virial speed changes which means that the number density, n_i, and relative speed, Δv_ik, changes for all particles and the full rate matrix should need an update; however, since n_i ∝ R(t)^-3 ∝ η^-3 and , by keeping track of the ηs it is not necessary to update all r_i(k)-elements every time step. Other quantities, like r (Eq. (18)), have simple scalings with η and η_K.

Finally we also need to update the cumulative distribution functions (CDFs) for picking the representative and physical particles. The CDFs for the physical particles with i ≠ j only need to be updated from element j and onwards since the elements r_i(k) with k ≠ j are unchanged. The CDF for the representative particle and the CDF for physical particle with i = j need a complete update. The result of this is that the time for updating scales as ~n × (n − j + 1).

Fig. 3 Collapse time as a function of solid planetesimal radius for clouds with initially monodisperse mm-, cm-, or dm-sized pebbles. Power-law fits to the data and the free-fall time of the initial cloud (tff ~ 25.1 yrs) are added to the plot. For small planetesimals (≲50 km) the collapse time is inversely proportional to the planetesimal radius and increases linearly with the pebble size. This is the same relation as we find in Appendix A where we assume bouncing as the only collisional outcome. For larger planetesimals (Rsolid ≳ 50 km) the collapse time drops rapidly as pebbles fragment and the number density and collision rates increases. The collapse time is, however, limited downwards by free-fall, so for even larger planetesimals with Rsolid ≳ 70 km the collapse time is equal to the free-fall time of the initial cloud.

2.4. Pooling of collisions

Because of high relative speeds some collisions in the simulations can lead to complete fragmentation. The collision rate between a large particle and monomers increases drastically because of the large number density of small particles (see Eq. (15), ). This, in turn, leads to very short time steps and many collisions with almost no impact on the system. We solve this issue by pooling up collisions if the mass ratio of the two particles is too small, f = m_k/m_i<f_crit. We pool collisions using the method suggested in Zsom & Dullemond (2008): we lower the probability for collisions with small mass ratios (projectile mass ≪target mass), but increase the number of physical collisions accordingly, .

The effect of the collision changes with pooling. The difference between 1 and X collisions is different for different collisional outcomes:

• Bouncing: everything as usual butδE → X·δE.

• Coagulation: same as for bouncing when it comes to the energy, but the mass and other physical properties of the representative particle changes as

• Fragmentation: for fragmentation we have two cases:

  • 1)

    Complete fragmentation of the representative particle with≤X collisions needed.

  • 2)

    Partial fragmentation even when adding up X collisions.

  In Case 1) we assume that Y = E_frag/E_coll collisions are required to completely fragment particle i and that the rest of the collisions result in bouncing with the fragment. The new representative particle will have monomer properties and the total energy loss will be δE = δE_frag + δE_bounce, where δE_frag is the number of physical particles in swarm i times the energy required to completely fragment one of them and δE_bounce is the energy lost from X − Y pooled collisions between a representative particle of monomer mass, m₀, and a physical particle of mass m_k (times the number of physical particles in the swarm, N_i). Case 2 is more likely and the outcome is calculated just like for one collision but where the collision energy has increased, E_coll → X·E_coll.

Fig. 4 Total energy, parametrized through η ≡ R(t) /R0, as a function of time for two clouds with, initially, cm-sized pebbles. Left panel: collapse of a cloud with Rsolid = 5 km; right panel: cloud with Rsolid = 50 km. The red lines show the actual size of the cloud, blue lines the equilibrium values (ηeq ≡ E0/E), and green lines the kinetic energy (ηK ≡ T0/T). The pink line shows the free-fall collapse of a cloud and the teal lines the analytic solution from Appendix A. We note that time is measured in units of the collapse time. Therefore the free-fall line is different for the two clouds. In the case of the low-mass cloud the collapse is slow and the only outcome of collisions is bouncing. Therefore it is in virial equilibrium until the very end of the collapse and follows the analytical solution (η = ηK = ηeq). The massive cloud, on the other hand, has initially higher collision speeds which lead to pebble fragmentation and rapid energy dissipation. At about η = 0.5 the energy dissipation becomes faster than free-fall, so the cloud does not get into virial equilibrium before the next pebble collision (the first term in Eq. (7)dominates over the second). This results in a transition to cold collapse (η is parallel to the free-fall curve) with subvirial speeds (ηK>ηeq).

3. Results

With our simulations we are interested in finding the evolution of the cloud properties during the collapse. We explicitly aim to measure the collapse time as a function of different cloud and pebble parameters, mainly the total mass of the cloud. We also follow the evolution of the size distribution of the pebbles. If the cloud is massive, then the relative speeds between the pebbles will be high and the collisions may result in fragmentation (see Fig. 1). From Eq. (1)we find that the virial speed in the initial state is proportional to the size of the planetesimal which will form, (22)Here we assumed that the initial cloud radius is equal to the Hill radius at its semi-major axis. The radius R_solid is the radius of a body with mass M and density equal to the monomer density. Eq. (22)implies that even if a massive cloud initially consists of cm-sized pebbles, these pebbles can be ground down to μm-sized dust during the collapse. For low-mass clouds the speeds are lower and collisions will only result in bouncing and the primordial pebbles survive the collapse.

3.1. Initial conditions and simulation setup

In our simulations we start out with a self-gravitating cloud of silicate pebbles with a monodisperse size distibution. The pebbles are assumed to be homogeneous spheres (density ρ_SiO2 ~ 2.5 g cm^-3) built up of μm-sized monomers. The reason for this density is that only the collisional outcome regions for silicates are known in details from laboratory experiments (Fig. 1). We assume that the initial radius of the cloud is equal to the Hill radius corresponding to the cloud mass. This means that the density of the cloud, and consequently the free-fall time, is independent of the cloud mass. To be able to use our zero-dimensional model we assume that the cloud is uniform, has no net angular momentum and always strives to get into virial equilibrium. We investigate the planetesimal formation at Kuiper belt distances from the Sun (at orbits with a semi-major axis equal to Pluto’s current orbit). Finally we also neglect the influence of any surrounding gas, as the friction timescale is much longer than the collision timescale in those dense clouds.

3.2. Collapsing pebble clouds

The range of cloud sizes we investigate corresponds to solid planetesimals with radii between 1 and 1000 km, from the smallest planetesimal sizes up to radii comparable to Pluto’s. In Fig. 3 we show the results of our simulations of the collapse time, the time at which the cloud reaches monomer density, as a function of planetesimal radius. In the figure we can see that for the simulations with cm-sized pebbles the collapse can be split into three regimes depending on the size of the planetesimal:

• 1) R_solid ≲ 50 km:

  The energy dissipation is dominated bybouncing collisions between the primordial pebbles. Collisionspeeds only reach values high enough for fragmentation to oc-cur in the very end of the collapse. The power-law fit shows us that which is the same dependence as we get from the analytic derivations when we assume that the particle size is constant and bouncing is the only collisional outcome (Appendix A). We find that the collapse time increases linearly with pebble size, which also agrees with the analytic model in Appendix A.

• 2) 50km ≲ R_solid ≲ 70 km:

  In this regime the cloud is massive enough (increased collision speeds) for collisions to result in pebble fragmentation early enough in the collapse to affect the collapse time. From Appendix A we have t_coll ∝ a/R_solid with constant particle size so the cloud collapses faster if the particles are smaller (increased collision rate). We find that the amount of fragmentation is larger and starts earlier (see Fig. 6) the more massive the cloud is (larger initial velocities). This results in a more rapid decrease of the collapse time with increasing cloud mass than in the previous regime.

• 3) R_solid ≳ 70 km:

  Collisions are so frequent (Eq. (15)) and dissipative (Eq. (9)) that the collapse is completely limited by the free-fall time of the cloud. Since the clouds initially have the same density, regardless of the cloud mass, above a solid radius of ~70 km all clouds have collapse times equal to the free-fall time.

For the simulations with mm-sized pebbles we only get the first and third regime. The reason for this is that the collapse gets limited by free-fall before it is massive enough for the relative speeds to be high enough for fragmenting collisions. We can also see in Fig. 3 that the collapse time of a cloud in regime 1 is a factor of 10 smaller with mm-sized pebbles than with cm-sized pebbles, as expected from the analytic derivations (t_coll ∝ a/R_solid).

Fig. 5 Virial and actual speeds from η and ηK as a function of time for three clouds of solid radii of 5 km (green and black), 50 km (blue and gold), and 500 km (cyan and red) with, initially, cm-sized pebbles. The low-mass cloud has a steadily increasing velocity because collisions are infrequent so that virialization can happen for each value of E. The massive cloud, on the other hand, dissipates energy too quickly and collapses cold with random pebble speeds much lower than the virial value. The initial virial speed is higher for the more massive planetesimal since the speed is a growing function of total mass (Eq. (1)).

The details of a collapse are different depending on the mass of the cloud. Figure 4 shows the collapse parameter, η, as a function of time for two different masses. The left plot follows the collapse of a cloud in regime 1 and the right plot a cloud in regime 2. If we would have plotted the collapse of a cloud in regime 3 it would have followed the free-fall curve closely. In the low-mass cloud we can see that the cloud collapses smoothly and the actual size (η, red curve) follows the equilibrium value (η_eq, blue curve). For a cloud of higher mass we run into the problem with virialization and at some point it cannot contract fast enough to keep up with the energy loss in the collisions. At this point η_eq<η and |δη| > |δη_max| (Eq. (5)) and the energy dissipation from this point on is so quick that the cloud free-falls. Hence the η-curve is parallel to the free-fall curve. Clouds of higher mass in regime 2 have their first fragmenting collision earlier in their collapse and therefore a shorter collapse time.

We find in the simulations that both the virialization and the energy dissipation in the collisions are important for the evolution of the kinetic energy. In Fig. 5 we plot the evolution of the speed of the particles which we get from Eq. (4), assuming that all particles have equal speeds. For the cloud in regime 1 (black curve) the speeds are steadily increasing thanks to the virialization of the cloud except for the very end of the collapse where virialization cannot keep up. For the cloud in regime 2 (yellow curve), on the other hand, the collapse is, as we have already seen, somewhat limited by free-fall and the speeds decrease thanks to the energy dissipation in the collisions (right panel in Fig. 4). In this case the first term in Eq. (7)dominates over the second term. This effect is greater the more massive the cloud is. The particles in the cloud in regime 3 (black curve) move with subvirial speeds very early in the collapse and the collision speeds reach values low enough for coagulation. We also see that the initial virial speed scales linearly with planetesimal radius because the clouds initially have a radius equal to their Hill radius and they are in virial equilibrium so v_vir,0 ∝ R₀ (Eq. (1)).

Fig. 6 Pebble mass fraction in the resulting planetesimals for simulations with, initially, cm-sized pebbles, as a function of the planetesimal radius, Rsolid. Particles are defined as pebbles if they have radii a> 1 mm. For more massive planetesimals the collision speeds are larger and the collisions can then result in fragmentation, decreasing the pebble mass fraction. However, massive clouds do not have time to virialize after each collision and the particles will move with subvirial speeds (see Fig. 5). This causes the number of fragmenting collisions to decrease as the collapse goes on, ensuring that a fraction of the primordial pebbles survive even in the most massive planetesimals.

Finally we are also interested in the sizes of the particles in the resulting planetesimal. If the planetesimal is built up by cm-sized pebbles it will be porous and have both low density and low internal strength. Such an object could be tidally broken up like the Shoemaker-Levy 9 comet (Asphaug & Benz 1996). If, on the other hand, the planetesimal is a mixture of pebbles and dust it can be packed very compactly and have a high density. In the case of low-mass clouds they collapse without any fragmentation and therefore become pebble piles with low density. This is in agreement with both the observed mass-density relation of Kuiper belt objects (low-mass Kuiper belt objects have low density, Brown 2013) and the model of comets as pebble piles with low internal strength (e.g. Skorov & Blum 2012; Blum et al. 2014). The density of planetesimals can, however, evolve as time goes on after the collapse through other physical mechanisms. Self-gravity can cause static compression and an increase in the density of the planetesimals (e.g. Kataoka et al. 2013). Another mechanism that could affect both the internal structure and the pebble mass fraction of the planetesimals is radioactive heating. Short-lived radioactive isotopes such as ²⁶Al could supply enough heat for differentiation to occur.

For a high-mass cloud fragmentation of pebbles is occurring during the collapse and the planetesimal can be more tightly packed and have a higher density, again consistent with the higher densities of larger Kuiper belt objects. In Fig. 6 we plot, for simulations with initially cm-sized pebbles, the pebble mass fraction (particles with radii a> 1 mm) in the resulting planetesimals as a function of the size of the planetesimals. First of all we see that planetesimals with radii R_solid ≲ 20 km are completely built up by pebbles. Planetesimals slightly larger, 20 km ≲ R_solid ≲ 50 km, are still in collapse regime 1 and follow the analytic approximations relatively well because the first fragmenting collision occurs late in the collapse and energy dissipation is dominated by bouncing collisions between cm-sized pebbles. Next we see that the pebble mass fraction decreases with increasing planetesimal size. However, even for the most massive planetesimals we still have a large fraction of the mass in pebbles. The reason that not all pebbles are ground down to dust is that the particles move with subvirial speeds in the end of the collapse (see Fig. 5) and the collision speeds are not high enough for fragmentation to occur. In the end of the simulations most collisions result in bouncing or coagulation.

The evolution of the particle size distribution in a simulation of a cloud collapsing into a 100 km-sized planetesimal is shown in Fig. 7 and can broadly be summarized as follows. The cloud starts with a monodisperse size distribution of cm-sized pebbles. As times goes on, pebbles collide with high enough velocities to fragment down to μm-sized dust. The subvirial collision speeds in the cold collapse phase ensures that a fraction of the pebbles survive. Collisions involving the dust can result in sticking and the dust coagulate and grow up to ~mm-sized pebbles again. In the end there is a roughly bimodal size distribution consisting of some primordial pebbles have survived the collapse and the growing dust that originally was ground down from the primordial pebbles.

Fig. 7 Particle size distribution in a collapsing pebble cloud with solid radius Rsolid = 100 km at four different times. The collision speeds are initially high enough for fragmentation to occur. Later in the collapse the particles move with subvirial speeds (see Fig. 5) and the amount of fragmentation decreases. The low collision speed in the cold collapse allow the dust to coagulate and grow to sizes of ~0.1 mm. One thing we also see in the figure is that, at the end of the collapse, a small fraction of the primordial pebbles have been subject to erosion through collisions with dust and decreased in size.

4. Conclusions and discussion

In this paper we have investigated the internal evolution of pebble clouds formed in protoplanetary disks by the streaming instability, although the results can be applied to other planetesimal formation models relying on particle concentration and gravitational collapse as well (e.g. Klahr & Bodenheimer 2003; Cuzzi et al. 2008; Dittrich et al. 2013). We have developed a numerical model to study the effect of dissipative collisions between pebbles leading to gravitational collapse of the cloud. The main results of our simulations can be summarized as follows.

The collapse times of the pebble clouds in our model are short. We assume that all clouds initially have a size equal to their Hill radius, giving a free-fall time of t_ff ~ 25.1 years at ~40 AU. The collapse time still varies depending on the cloud properties. Massive clouds (R_solid ≳ 70 km) collapse, roughly, on the free-fall time of the cloud and even for low-mass clouds (R_solid ~ km) the collapse time is only a few tens of orbits. The size of the pebbles affects the collapse time, with more rapid collapse for smaller pebbles (because of their larger collision surface). These collapse timescales can be compared to the formation times of the pebble clouds which are just a few orbits (Johansen et al. 2009).

The clouds do not get into virial equilibrium immediately after a dissipative collision but require some time to fall to their new (smaller) desired size. The result is that the particles in high-mass clouds will move at subvirial speeds, reducing their collision speeds significantly and preventing widespread fragmentation. For a low-mass cloud this does not occur during most of the collapse, since the collision rates are low and the cloud has time to virialize between two collisions. It is only at the end of the collapse, when the density of the cloud is very high, that the particles in low-mass clouds dissipate energy quickly enough that they will move with subvirial speeds. A high-mass cloud, on the other hand, enters the cold collapse phase already in the beginning of the contraction. The result is that massive clouds collapse coldly and after an initial burst of pebble fragmentation the number of fragmenting collisions decreases and the dust can start to grow to reform larger and larger aggregates (up to ~mm).

In the collapse of low-mass clouds (R_solid ≲ 20 km) the collision speeds never reach the fragmentation limit and they will be completely made up of the primordial pebbles. In more massive clouds (R_solid ≳ 20 km) the collision speeds are, at some point in the collapse, large enough for fragmentation to occur in the pebble collisions. In all planetesimals some of the primordial pebbles survive. For the mid-size planetesimals (20km ≲ R_solid ≲ 50km) the first fragmenting collision occurs late in the collapse and not all pebbles have time to be ground down. In the most massive clouds even the first collisions are fragmenting but slow virialization of the cloud causes the collision speeds later in the collapse to be smaller and a significant fraction of the pebbles survive. This formation process causes an observable difference between a low- and a high-mass planetesimal, namely the bulk density. Low-mass planetesimals can be thought of as pebble piles and they will be very loosely packed and have a low density. This will also make them easily tidally disrupted observed in Jupiter’s tidal disruption of the comet Shoemaker-Levy 9. More massive planetesimals, containing both dust and pebbles, are able to become very tightly packed and therefore achieve a high density. This relationship has been suggested for objects in the Kuiper belt (Brown 2013) but we stress that there are more parameters that affect the density, such as composition, differentiation through radioactive heating and collisional impacts.

We use the results from laboratory experiments of collisions between silicate aggregates even if Kuiper belt objects are mainly composed of ice. One difference is that ice particles could be very fluffy and survive collisions at higher speeds compared to silicates (e.g. Wada et al. 2009, 2013). However, even in our most massive clouds the collision speeds do not reach values over ~10 m s^-1 and the speed decreases rapidly as the collapse progresses (red line in Fig. 5). The collapse times, on the other hand, could very well be affected by our choice of material. In Appendix A we can see that for a given cloud mass and pebble size the collapse time decreases for ice since the material density is lower (the collision rates increases) but the main conclusion that pebbles survive the collapse would not change.

The next step in our investigation will be to increase the dimensionality of the simulations. Our current numerical model is zero-dimensional in the sense that we only keep track of the cloud size and collisions are considered via a Monte Carlo model. The 0D approach is not completely physically correct since the pebble cloud should initially have some rotation which will prevent the direct collapse. However, we expect that the main findings of the paper will be true also in a three-dimensional model with rotation, since gravitational instabilities in the rotating pebble disk should lead to angular momentum transport and the excitation of relative speeds which can drive a collapse similar to the zero-dimensional model. Hydrodynamical two- and three-dimensional models will be presented in a subsequent paper.

Acknowledgments

K.W.J. and A.J. were supported by the European Research Council under ERC Starting Grant agreement 278675-PEBBLE2PLANET. A.J. was also supported by the Swedish Research Council (grant 2010-3710) and by the Knut and Alice Wallenberg Foundation. K.W.J. and A.J. acknowledges the support from the Royal Physiographic Society in Lund for grants to purchase computer hardware to run the simulations on. We are grateful for stimulating discussions with Jürgen Blum, Hubert Klahr, Tristen Hayfield, Patrick Glaschke, Thomas Henning, and Andrew Youdin. We also thank the referee Hiroshi Kobayashi for many insightful comments that helped improve the manuscript.

References

Appendix A: Analytic derivation of the collapse of a particle cloud

As mentioned in the main text, it is possible to derive the evolution of a collapsing cloud analytically under some assumptions. First of all we assume a uniform density and no net rotation so that we can make the model zero-dimensional. Next we assume that we have equal particle sizes, equal speeds (), and that the only outcome of a collision is bouncing. The last point is only valid for low-mass clouds which exhibit low collision speeds. A final assumption is that the cloud virializes immediately after a collision. This, we will later see, makes the relative speeds too high and the collapse time too small for clouds that are not in regime 1 (see Sect. 3.2).

We estimate the evolution of the energy of the cloud with the equation (A.1)where δE is the change in energy per collision and r_coll is the collision rate (number of collisions per unit time). The energy change per collision is (from Eq. (9)) (A.2)where m is the mass of each individual particle. The collisions are not necessarily head-on which makes the effective collision speed smaller and reduces the average energy dissipated with a factor of 1 / 2. To get the collision rate of the system we look at the rate with which two equal-sized particles in a volume V collide: (A.3)where σ_ik is the geometrical cross-section of particle i and k, a is the particle radius, and R is the radius of the cloud. The total rate is the sum of the rates of all possible collisions (a particle cannot collide with itself and we cannot count both r_ik and r_ki as it is the same collision), and setting r_ik constant for all collision pairs, we get (A.4)The last step is valid since the number of particles N_p ≫ 1. We then write the rate as (A.5)Here M is the total mass of the cloud and R_solid is the radius of a planetesimal with mass M and density equal to the internal density of the particles. Now Eq. (A.1)becomes (A.6)To get Δv and R as functions of the total energy we use Eq. (1)where the change from to Δv comes from the fact that the relative speed is on average times larger than the individual speeds. Inserting Eqs. (A.7), (A.8)into Eq. (A.6)yields (A.9)Here K = 3.5K′ is a measure of the rate with which the energy changes and Q = 3.5Q′ is a measure of the initial energy. To get the value of the parameter Q we look at the initial energy, E₀ = E(t = 0), and using Eq. (1)we obtain (A.10)where R₀ is the initial radius of the cloud. A quick look at Eq. (A.9)shows us that if t = t_coll ≡ Q/K then − E → ∞, i.e. a complete collapse (R → 0 and ). Finally we assume that the initial radius of the cloud is equal to the Hill radius at the semi-major axis of the cloud (A.11)where d is the semi-major axis of the cloud’s orbit. Now we rewrite Eq. (A.9)as (A.12)By using Eqs. (A.2), (A.5), (A.7)−(A.11), assuming silicate pebbles, and placing the cloud at Pluto’s distance from the Sun we can find an expression for the collapse time (A.13)

All Figures

	Fig. 1 Collisional outcomes as a function of projectile radius, a, collision speed, Δv, and relative particle size, f ≡ m2/m1. Green regions indicate coagulation, yellow regions bouncing, and red regions fragmentation. Outcome regions taken from Güttler et al. (2010).
In the text

Fig. 2 In our model pebbles are composed of μm-sized monomer particles. In collisions resulting in partial fragmentation monomers are removed one by one to simulate erosion. The binding energy of the fragment, , is larger (closer to zero) than the binding energy of the original particle, Ebind, since the fragment consists of fewer monomers. The number of monomers released, Nmon, is equal the number of contact surfaces the collision energy can break divided by three contact surfaces per monomer, Nmon = Ecoll/(3Eroll).

In the text

Fig. 3 Collapse time as a function of solid planetesimal radius for clouds with initially monodisperse mm-, cm-, or dm-sized pebbles. Power-law fits to the data and the free-fall time of the initial cloud (tff ~ 25.1 yrs) are added to the plot. For small planetesimals (≲50 km) the collapse time is inversely proportional to the planetesimal radius and increases linearly with the pebble size. This is the same relation as we find in Appendix A where we assume bouncing as the only collisional outcome. For larger planetesimals (Rsolid ≳ 50 km) the collapse time drops rapidly as pebbles fragment and the number density and collision rates increases. The collapse time is, however, limited downwards by free-fall, so for even larger planetesimals with Rsolid ≳ 70 km the collapse time is equal to the free-fall time of the initial cloud.

In the text

Fig. 4 Total energy, parametrized through η ≡ R(t) /R0, as a function of time for two clouds with, initially, cm-sized pebbles. Left panel: collapse of a cloud with Rsolid = 5 km; right panel: cloud with Rsolid = 50 km. The red lines show the actual size of the cloud, blue lines the equilibrium values (ηeq ≡ E0/E), and green lines the kinetic energy (ηK ≡ T0/T). The pink line shows the free-fall collapse of a cloud and the teal lines the analytic solution from Appendix A. We note that time is measured in units of the collapse time. Therefore the free-fall line is different for the two clouds. In the case of the low-mass cloud the collapse is slow and the only outcome of collisions is bouncing. Therefore it is in virial equilibrium until the very end of the collapse and follows the analytical solution (η = ηK = ηeq). The massive cloud, on the other hand, has initially higher collision speeds which lead to pebble fragmentation and rapid energy dissipation. At about η = 0.5 the energy dissipation becomes faster than free-fall, so the cloud does not get into virial equilibrium before the next pebble collision (the first term in Eq. (7)dominates over the second). This results in a transition to cold collapse (η is parallel to the free-fall curve) with subvirial speeds (ηK>ηeq).

In the text

Fig. 5 Virial and actual speeds from η and ηK as a function of time for three clouds of solid radii of 5 km (green and black), 50 km (blue and gold), and 500 km (cyan and red) with, initially, cm-sized pebbles. The low-mass cloud has a steadily increasing velocity because collisions are infrequent so that virialization can happen for each value of E. The massive cloud, on the other hand, dissipates energy too quickly and collapses cold with random pebble speeds much lower than the virial value. The initial virial speed is higher for the more massive planetesimal since the speed is a growing function of total mass (Eq. (1)).

In the text

Fig. 6 Pebble mass fraction in the resulting planetesimals for simulations with, initially, cm-sized pebbles, as a function of the planetesimal radius, Rsolid. Particles are defined as pebbles if they have radii a> 1 mm. For more massive planetesimals the collision speeds are larger and the collisions can then result in fragmentation, decreasing the pebble mass fraction. However, massive clouds do not have time to virialize after each collision and the particles will move with subvirial speeds (see Fig. 5). This causes the number of fragmenting collisions to decrease as the collapse goes on, ensuring that a fraction of the primordial pebbles survive even in the most massive planetesimals.

In the text

Fig. 7 Particle size distribution in a collapsing pebble cloud with solid radius Rsolid = 100 km at four different times. The collision speeds are initially high enough for fragmentation to occur. Later in the collapse the particles move with subvirial speeds (see Fig. 5) and the amount of fragmentation decreases. The low collision speed in the cold collapse allow the dust to coagulate and grow to sizes of ~0.1 mm. One thing we also see in the figure is that, at the end of the collapse, a small fraction of the primordial pebbles have been subject to erosion through collisions with dust and decreased in size.

In the text