Issue 
A&A
Volume 570, October 2014



Article Number  A47  
Number of page(s)  10  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201424369  
Published online  15 October 2014 
Formation of pebblepile planetesimals
Lund ObservatoryDepartment of Astronomy and Theoretical Physics, Lund
University,
Box 43,
221 00
Lund,
Sweden
email:
kalle@astro.lu.se
Received:
11
June
2014
Accepted:
1
September
2014
Asteroids and Kuiper belt objects are remnant planetesimals from the epoch of planet formation. The first stage of planet formation is the accumulation of dust and ice grains into mm and cmsized pebbles. These pebbles can clump together through the streaming instability and form gravitationally bound pebble clouds. Pebbles inside such a cloud will undergo mutual collisions, dissipating energy into heat. As the cloud loses energy, it gradually contracts towards solid density. We model this process and investigate two important properties of the collapse: (i) the collapse timescale and (ii) the temporal evolution of the pebble size distribution. Our numerical model of the pebble cloud is zerodimensional and treats collisions with a statistical method. We find that planetesimals with radii larger than ~100 km collapse on the freefall timescale of ~25 years. Lowermass clouds have longer pebble collision timescales and collapse much more slowly, with collapse times of a few hundred years for 10 km scale planetesimals and a few thousand years for 1 km scale planetesimals. The mass of the pebble cloud also determines the interior structure of the resulting planetesimal. The pebble collision speeds in lowmass clouds are below the threshold for fragmentation, forming pebblepile planetesimals consisting of the primordial pebbles from the protoplanetary disk. Planetesimals above 100 km in radius, on the other hand, consist of mixtures of dust (pebble fragments) and pebbles which have undergone substantial collisions with dust and other pebbles. The Rosetta mission to the comet 67P/ChuryumovGerasimenko and the New Horizons mission to Pluto will provide valuable information about the structure of planetesimals in the solar system. Our model predicts that 67P is a pebblepile planetesimal consisting of primordial pebbles from the solar nebula, while the pebbles in the cloud which contracted to form Pluto must have been ground down substantially during the collapse.
Key words: methods: analytical / methods: numerical / minor planets, asteroids: general / planets and satellites: formation / comets: general
© ESO, 2014
1. Introduction
In the process of forming planets the initially μmsized dust and ice particles in the primordial protoplanetary disk grow to planets of sizes up to ~10^{4} 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 KelvinHelmholz turbulence in the midplane 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 coformation explanation is not possible in models of binary formation through threebody encounters (Goldreich et al. 2002) in the hierarchical coagulation formation mechanism.
Nesvorný et al. (2010) modelled the particle cloud with an Nbody 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. Lowmass clouds experience no fragmentation and the planetesimal will be formed from primordial pebbles. The pebbles in highmass 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 lowmass pebble cloud where the collision speeds never reach the fragmentation limit. Studies of the tidal breakup of the comet ShoemakerLevy 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/ChuryumovGerasimenko (Bauer et al. 2012) show that mmsized pebbles are released from the surface of the comet. Skorov & Blum (2012) model this pebble release process with a twolayered comet nucleus consisting of an icefree 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 cmsized pebbles. The interior structure of comets will also be investigated by the Rosetta mission which will put a lander on 67P/ChuryumovGerasimenko 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 PlutoCharon system in July 2015 (Stern 2008). After a flyby 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 flybys. As no currently known Kuiper belt objects are within the range of New Horizons after the flyby 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 zerodimensional 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 lowmass 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 zerodimensional, 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_{0}, is (1)Here T_{0} 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_{0} is the initial potential energy, G is the gravitational constant, and R_{0} 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 freefall 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 freefalling from R_{0}. This gives (5)where t_{ff,0} is the initial freefall time. Here we have used the freefall speed at R starting from R_{0}, (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 ≡ m_{2}/m_{1}. 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_{2}/m_{1}, of the particles. The collisions are split into two groups:

1.
two similarsized 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:
Δv (m s^{1})  f ≥ 0.1  f< 0.1 
Δv<v_{stick}  C  C 
v_{stick} ≤ Δv< 1  B  B 
1 ≤ Δv< 25  F  C 
Δv ≥ 25  F  F 
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_{1}m_{2}/ (m_{1} + m_{2}) is the reduced mass, a_{0} is the monomer radius, and F_{roll} is the rolling force of the monomers^{1}. 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 nonzero 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 μmsized 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, E_{bind}, since the fragment consists of fewer monomers. The number of monomers released, N_{mon}, is equal the number of contact surfaces the collision energy can break divided by three contact surfaces per monomer, N_{mon} = E_{coll}/^{(}3E_{roll}^{)}. 
The effect of fragmentation is a bit different. In a collision between two particles with mass m_{1} and m_{2} we have, in the centreofmass frame, an available collision energy of (10)Next we assume that all particles are composed of μmsized monomers of mass m_{0}. 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 μmsized 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_{1} and a_{2}) of the two particles with higher probability for larger impact parameters (larger area), (14)The effect of a nonzero 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 cmsized 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 Plutomass planetesimal split into cmsized pebbles results in N ~ 10^{24} 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 crosssection 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 onedimensional 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_{0}. If, on the other hand, we have partial fragmentation, we pick the new mass from a massweighted 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_{0}. 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 dmsized pebbles. Powerlaw fits to the data and the freefall time of the initial cloud (t_{ff} ~ 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 (R_{solid} ≳ 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 freefall, so for even larger planetesimals with R_{solid} ≳ 70 km the collapse time is equal to the freefall 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_{0}, 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}.

1)
Fig. 4
Total energy, parametrized through η ≡ R(t) /R_{0}, as a function of time for two clouds with, initially, cmsized pebbles. Left panel: collapse of a cloud with R_{solid} = 5 km; right panel: cloud with R_{solid} = 50 km. The red lines show the actual size of the cloud, blue lines the equilibrium values (η_{eq} ≡ E_{0}/E), and green lines the kinetic energy (η_{K} ≡ T_{0}/T). The pink line shows the freefall 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 freefall line is different for the two clouds. In the case of the lowmass 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 freefall, 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 freefall 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 semimajor 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 cmsized pebbles, these pebbles can be ground down to μmsized dust during the collapse. For lowmass 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 selfgravitating cloud of silicate pebbles with a monodisperse size distibution. The pebbles are assumed to be homogeneous spheres (density ρ_{SiO}_{2} ~ 2.5 g cm^{3}) built up of μmsized 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 freefall time, is independent of the cloud mass. To be able to use our zerodimensional 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 semimajor 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 cmsized 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 occur in the very end of the collapse. The powerlaw 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 freefall 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 freefall time.
For the simulations with mmsized pebbles we only get the first and third regime. The reason for this is that the collapse gets limited by freefall 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 mmsized pebbles than with cmsized 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, cmsized pebbles. The lowmass 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 freefall curve closely. In the lowmass 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 freefalls. Hence the ηcurve is parallel to the freefall 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 freefall 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_{0} (Eq. (1)).
Fig. 6
Pebble mass fraction in the resulting planetesimals for simulations with, initially, cmsized pebbles, as a function of the planetesimal radius, R_{solid}. 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 cmsized pebbles it will be porous and have both low density and low internal strength. Such an object could be tidally broken up like the ShoemakerLevy 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 lowmass clouds they collapse without any fragmentation and therefore become pebble piles with low density. This is in agreement with both the observed massdensity relation of Kuiper belt objects (lowmass 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. Selfgravity 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. Shortlived radioactive isotopes such as ^{26}Al could supply enough heat for differentiation to occur.
For a highmass 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 cmsized 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 cmsized 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 kmsized planetesimal is shown in Fig. 7 and can broadly be summarized as follows. The cloud starts with a monodisperse size distribution of cmsized pebbles. As times goes on, pebbles collide with high enough velocities to fragment down to μmsized 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 ~mmsized 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 R_{solid} = 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 freefall 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 freefall time of the cloud and even for lowmass 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 highmass clouds will move at subvirial speeds, reducing their collision speeds significantly and preventing widespread fragmentation. For a lowmass 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 lowmass clouds dissipate energy quickly enough that they will move with subvirial speeds. A highmass 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 lowmass 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 midsize 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 highmass planetesimal, namely the bulk density. Lowmass 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 ShoemakerLevy 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 zerodimensional 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 threedimensional 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 zerodimensional model. Hydrodynamical two and threedimensional models will be presented in a subsequent paper.
We use F_{roll} = (8.5 ± 1.6) × 10^{10} N for SiO_{2} spheres with a_{0} ~ 1μm from Heim et al. (1999).
Acknowledgments
K.W.J. and A.J. were supported by the European Research Council under ERC Starting Grant agreement 278675PEBBLE2PLANET. A.J. was also supported by the Swedish Research Council (grant 20103710) 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
 Asphaug, E., & Benz, W. 1996, Icarus, 121, 225 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Barrat, A., & Trizac, E. 2002, Granular Matter 4, 57 [CrossRef] [Google Scholar]
 Bauer, J. M., Kramer, E., Mainzer, A. K., et al. 2012, ApJ, 758, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Benecchi, S. D., Noll, K. S., Grundy, W. M., et al. 2009, Icarus, 200, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2014, Icarus, 235, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, M. E. 2013, ApJ, 778, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Buie, M. W., Spencer, J. R., Parker, A. H., et al. 2012, LPI Contributions, 1667, 6430 [NASA ADS] [Google Scholar]
 Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., Gutiérrez, P. J., & Rickman, H. 2007, Icarus, 187, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Dittrich, K., Klahr, H., & Johansen, A. 2013, ApJ, 763, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., Lithwick, Y., & Sari, R. 2002, Nature, 420, 643 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heim, L.O., Blum, J., Preuss, M., & Butt, H.J. 1999, Phys. Rev. Lett., 83, 3328 [NASA ADS] [CrossRef] [Google Scholar]
 Higa, M., Arakawa, M., & Maeno, N. 1996, Planet. Space Sci., 44, 917 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low, M.M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, & Th. Henning, submitted [arXiv:1402.1344] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 554, A4 [Google Scholar]
 Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Safronov, V. S. 1969, Evoliutsiia Doplanetnogo Oblaka (English transl.: Evolution of the Protoplanetary Cloud and Formation of Earth and the Planets, NASA Tech. Transl. F677, Jerusalem: Israel Sci. Transl., 1972) [Google Scholar]
 Skorov, Y., & Blum, J. 2012, Icarus, 221, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Stern, S. A. 2008, Space Sci. Rev., 140, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Weissman, P. R., Asphaug, E., & Lowry, S. C. 2004, Comets II (University of Arizona Press), 337 [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
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 zerodimensional. 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 lowmass 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 headon 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 equalsized particles in a volume V collide: (A.3)where σ_{ik} is the geometrical crosssection 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_{0} = E(t = 0), and using Eq. (1)we obtain (A.10)where R_{0} 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 semimajor axis of the cloud (A.11)where d is the semimajor 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 ≡ m_{2}/m_{1}. 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 μmsized 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, E_{bind}, since the fragment consists of fewer monomers. The number of monomers released, N_{mon}, is equal the number of contact surfaces the collision energy can break divided by three contact surfaces per monomer, N_{mon} = E_{coll}/^{(}3E_{roll}^{)}. 

In the text 
Fig. 3
Collapse time as a function of solid planetesimal radius for clouds with initially monodisperse mm, cm, or dmsized pebbles. Powerlaw fits to the data and the freefall time of the initial cloud (t_{ff} ~ 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 (R_{solid} ≳ 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 freefall, so for even larger planetesimals with R_{solid} ≳ 70 km the collapse time is equal to the freefall time of the initial cloud. 

In the text 
Fig. 4
Total energy, parametrized through η ≡ R(t) /R_{0}, as a function of time for two clouds with, initially, cmsized pebbles. Left panel: collapse of a cloud with R_{solid} = 5 km; right panel: cloud with R_{solid} = 50 km. The red lines show the actual size of the cloud, blue lines the equilibrium values (η_{eq} ≡ E_{0}/E), and green lines the kinetic energy (η_{K} ≡ T_{0}/T). The pink line shows the freefall 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 freefall line is different for the two clouds. In the case of the lowmass 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 freefall, 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 freefall 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, cmsized pebbles. The lowmass 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, cmsized pebbles, as a function of the planetesimal radius, R_{solid}. 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 R_{solid} = 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.