Radially resolved simulations of collapsing pebble clouds in protoplanetary discs

TL;DR

This study models the collapse of pebble clouds in protoplanetary discs to understand the internal structure of comet-sized planetesimals, highlighting the influence of gas drag and collision outcomes on porosity and pebble size distribution.

Contribution

It introduces a model allowing differential contraction of pebble sub-clouds and incorporates gas drag effects, advancing understanding of planetesimal internal structure formation.

Findings

Comets are porous pebble-piles with depth-dependent pebble sizes.

Surface layers contain a mix of primordial pebbles and fragments.

Interior regions consist mainly of primordial pebbles with higher porosity.

Abstract

We study the collapse of pebble clouds with a statistical model to find the internal structure of comet-sized planetesimals. Pebble-pebble collisions occur during the collapse and the outcome of these collisions affect the resulting structure of the planetesimal. We expand our previous models by allowing the individual pebble sub-clouds to contract at different rates and by including the effect of gas drag on the contraction speed and in energy dissipation. Our results yield comets that are porous pebble-piles with particle sizes varying with depth. In the surface layers there is a mixture of primordial pebbles and pebble fragments. The interior, on the other hand, consists only of primordial pebbles with a narrower size distribution, yielding higher porosity there. Our results imply that the gas in the protoplanetary disc plays an important role in determining the radial distribution of pebble sizes and porosity inside planetesimals.

Equations

$${\dot{M}}_{\text{c}}=\sum _i\pi R_{\text{c}}^2{v}_i{n}_i{m}_i,$$

$$r_{ik}$$

$$r_{ki}$$

$$\delta E_{ik}$$

$$\delta t$$

$$E_{i,\text{tot}}$$

$$v_{i,\text{ff}}=\sqrt{\frac{2G{M}_{i,\text{enc}}}{R_{i,0}}}\sqrt{\frac{R_{i,0}-R_i}{R_i}},$$

$$v_{i,\text{t}}={\tau }_{i,\text{f}}{F}_{i,\text{g}}=\frac{\pi {\rho }_{\text{s}}{a}_i}{2{\Sigma }_{\text{d}}{\Omega }_{\text{K}}}\frac{G{M}_{i,\text{enc}}}{R_i^2},$$

$${\dot{v}}_{i,\text{d}}=-\frac{v_i}{{\tau }_{i,\text{f}}}.$$

$$\text{d}N(a)$$

$$M_{i,\text{enc}}$$

$$m_{ik,\text{enc}}$$

$$U_i$$

$$U_{ik}$$

$${\Phi }_{ik}(r)$$

$$⟹U_{ik}$$

$$2E_{i,\text{tot}}$$

$$U_{i,\text{eq}}$$

$$-\frac{3G{M}_{\text{c}}{m}_i}{2R_{i,\text{eq}}}$$

$$=\underset{{\Sigma }_1}{\underbrace{\frac{3G{m}_i}{10}\sum _{k\ne i}\frac{m_k}{R_k^3}}}{R}_{i,\text{eq}}^2-\underset{{\Sigma }_2}{\underbrace{3G{m}_i\left(\frac{M_{\text{c}}}{2}+\frac{2m_i}{5}\right)}}\frac{1}{R_{i,\text{eq}}}$$

$$-\underset{{\Sigma }_3}{\underbrace{\frac{3G{m}_i}{2}\sum _{k\ne i}\frac{m_k}{R_k}}},$$

$$⟹0$$

$$U_{i,\text{eq}}$$

$$-\frac{3G{M}_{\text{c}}{m}_i}{2R_{i,\text{eq}}}$$

$$=-\underset{{\Sigma }_4}{\underbrace{3G{m}_i\left(\sum _{k\ne i}\frac{m_k}{2}+\frac{2m_i}{5}+\frac{M_{\text{c}}}{2}\right)}}\frac{1}{R_{i,\text{eq}}}$$

$$+\underset{{\Sigma }_5}{\underbrace{\frac{3G{m}_i}{10}\sum _{k\ne i}{m}_k{R}_k^2}}\cdot \frac{1}{R_{i,\text{eq}}^3},$$

$$⟹0$$

Full text

††pubyear: 2017††pagerange: Radially resolved simulations of collapsing pebble clouds in protoplanetary discs–A

1 Introduction

Planet formation occurs in protoplanetary discs surrounding newborn stars by the growth of, initially, $\mu$m-sized solid dust and ice particles up to planets with characteristic sizes $10^{4}-10^{5}$ km (Safronov, 1969). Through collisions and coagulation by contact forces small dust particles can grow to larger and larger sizes (Blum & Wurm, 2008). At some size ($\sim$mm-cm) the sticking capabilities of the particles become too poor for coagulation and collisions rather result in bouncing and compactification. Instead, growth can continue by concentration of pebbles in turbulent gas and the formation of gravitationally bound pebble clouds (see review by Johansen et al., 2014). Such clouds can later experience a gravitational collapse and form solid planetesimals with sizes of order $\sim$10-100 km (Johansen et al., 2015; Simon et al., 2016).

The ESA space probe Rosetta has, with its visit to and observations of the Jupiter family comet 67P/Churyumov–Gerasimenko (hereafter 67P), taught us a lot about the structure and interior of comets. Observations of the surface of 67P find a ‘goosebump’ structure inside deep pits, which has been suggested to be primordial m-sized pebbles that formed the comet (Sierks et al., 2015). Images from the Philae lander find a pebble size scale closer to 1 cm (Mottola et al., 2015) more in agreement with particle coagulation limits (Birnstiel et al., 2012). Observations of shape and gravity field result in a bulk density of $\sim$530 kg m*-3*. Radar tomography measurements of the comet core (Pätzold et al., 2016; Kofman et al., 2015) suggest that 67P is homogeneous on length scales $<$3 m and very porous (70-75% depending on assumed dust-to-ice ratio). These measurements overall support the picture that 67P is a porous pile of pebbles from the solar protoplanetary disc.

The models of collapsing pebble clouds used in our previous work (Wahlberg Jansson et al., 2017, hereafter WJ+17) assume that pebble clouds are homogeneous, with uniform internal density, and non-rotating. The collapse itself is driven by loss of kinetic energy in inelastic collisions between pebbles. Loss of energy leads to contraction of the cloud and higher particle speeds when it tries to return to virial equilibrium (Sec. 2). These models have resulting planetesimals with uniform particle size distribution throughout the interior, from the centre to the surface. A key result of these models is that the collapse time-scale is a function of both cloud mass and initial particle sizes (Wahlberg Jansson & Johansen, 2014, hereafter WJJ14). Because of their larger total surface area, small particles dissipate energy more efficiently than large particles, resulting in a faster collapse. The cloud mass also determines the magnitudes of pebble collision speeds. For low-mass pebble clouds, speeds never reach high enough values for fragmenting collisions to occur and the cloud collapses into a porous pebble-pile This would be the case for planetesimals up to $\sim$20 km in radius, covering the size range for most comet cores.

However, in previous models we assumed, for simplicity, that protoplanetary disc gas can be ignored during the collapse. Particles of different sizes are affected differently when travelling through gas, so the gas should influence the details of the cloud collapse process. In this paper we extend the previous calculations by adding the effect of gas on the collapsing pebble cloud. We also expand the collapse model to allow the individual particle clouds to have different physical sizes. This allows us to calculate the interior structure of the resulting planetesimal, particularly how the pebble size distribution varies with distance from the planetesimal centre. Bouncing collisions during the cloud collapse cause porous pebbles to become more compact. This leads to resulting pebble sizes to not be the same as the initial, even without fragmenting collisions Lorek et al. (2016). This evolution of particle density is, however, not taken into account in our simulations.

The paper is organized as follows. The radially resolved pebble cloud model is described in Section 2. In Section 2.4 our implementation of the effect of disc gas is summarized. Our simulations of collapsing pebble clouds are presented in Section 3 and the results are discussed in Section 4. In Appendix A we present an analytic description of the calculation of potential energies in a pebble cloud undergoing a gravitational collapse.

2 Model

We model the formation of a planetesimal through the collapse of a gravitationally bound pebble cloud. Such a cloud can form inside a protoplanetary disc through e.g. the streaming instability (Youdin & Goodman, 2005; Johansen et al., 2009; Bai & Stone, 2010). In this cloud inelastic pebble-pebble collisions occur and dissipate kinetic energy from the cloud. This causes the cloud to contract and release kinetic energy to the pebbles due to the negative heat capacity property of self-gravitating systems. Increased pebble speeds and number density leads to more frequent collisions, higher energy dissipation rate and eventually a runaway collapse, the gravothermal catastrophe. Here we present an expansion of our previous numerical model (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17) to also include an investigation of the effect of gas on the collapse process and the radial structure inside the resulting planetesimals.

2.1 Model of a radially resolved gravitationally bound pebble cloud

In previous work (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17) we have assumed that all particles, regardless of size, follow each other as the cloud collapses. We keep track of the kinetic and potential energy of the entire cloud to find the radius, density, collision speeds etc. for the cloud and particle collisions. Effectively this approach spreads all the physical particles homogeneously over the entire cloud.

The expanded model includes the presence of the gas and the assignment of individual sizes to each collapsing sub-cloud of pebbles. The cloud is still modelled as a spherical, non-rotating cloud of pebbles but compared to the previous pebble cloud model we keep track of the radial distribution of pebbles and also allow for an initial pebble size distribution. The cloud is modelled as a number, $N_{\text{{p}}}$, of overlapping, spherical sub-clouds that together form a pebble cloud collapsing into a solid planetesimal (see Fig. 1). All sub-clouds are assumed to, at all times, seek virial equilibrium.

We assume that each representative particle (Zsom & Dullemond, 2008) is its own collapsing sub-cloud and has its own kinetic energy and size to describe its properties. Using multiple sub-clouds means that representative particles are affected by the contraction of another particle swarm. The enclosed mass of one sub-cloud increases if another, larger, sub-cloud contracts (increased density). This results in a change in potential energy. We take care of this by calculating the total potential energy of each cloud and update the size of all clouds at every time step.

The collapsing cloud will eventually reach the solid density in the centre. We replace the first particle to reach solid density with a sink particle at rest in the centre of the cloud. The sink particle (core) then accretes particles from the other swarms. The core mass, $M_{\text{{c}}}$, increases continuously as

[TABLE]

where $R_{\text{{c}}}$ is the core radius, $v_{i}$ is the random particle speed in swarm $i$, $n_{i}$ is the particle number density and $m_{i}$ is the mass of those particles. This way we keep track of the particle size distribution at different depths inside the final planetesimal. Individual particle speeds are independent on the physical size of the core ($R_{\text{{c}}}$) and only at the very end of the collapse, when most of the planetesimal mass has been accreted, is the escape speed from the planetesimal larger than the virial speed of the cloud. Therefore we neglect gravitational focusing in Eq. (1).

2.2 Particle collisions

In Fig. 1 we show the collision rates for two swarms of particles with different cloud sizes ($R_{i}<R_{k}$). The blue cloud is the denser $i$-swarm and the red cloud is the less dense $k$-swarm. The collision rate for $i$ being the representative particle is unchanged from the old model because $k$ encloses $i$ and we are only interested in the density of the $k$-swarm. The collision rate for $k$ being representative particle, on the other hand, has changed. Now we can only consider the red particles within $i$’s radius. This means that a factor $\left(R_{i}/R_{k}\right)^{3}$ has been added (the volume fraction of the small cloud in the large),

[TABLE]

This way we can continue as before and calculate all collision rates, calculate the time until next collision, $\delta t$, and which particles collide [Eqs. (17)-(21) in WJJ14]. The outcome (coagulation, bouncing, mass transfer and fragmentation) and energy lost, $\delta E_{ik}$, in a collision as function of relative speed and particle sizes are calculated as in WJ+17 with the results from experimental studies (Güttler et al., 2010; Bukhari Syed et al., 2017). In a collision a representative particle and hence a single sub-cloud can dissipate most of its energy in a single collision (a head-on collision with low coefficient of restitution). This is unrealistic but caused no problems in the homogeneous simulations in WJJ14 and WJ+17. In a collsion, a small fraction ($1/N_{\text{{p}}}$) of the particles lost kinetic energy. At the end of the time step, however, all particles were artificially given the same kinetic energy. In the new model, however, each sub-cloud keeps its energy for itself and losing all, or most, of its kinetic energy will cause free-fall of the sub-cloud after just one collision. Such a collision is unlikely but will happen when a large number of collisions occur. To resolve this ‘coarse-grained’ behaviour we add a ‘smoothing parameter’, $\zeta\in\left[0,1\right]$, which decreases the amount of energy dissipation per time step. The assumption, similar to the old model but implemented now for each individual sub-cloud, is that only a fraction of the particles (in the sub-cloud) collides and dissipate energy but all those particles share the same energy. This smoothing corresponds to decreasing the number of physical collisions that actually occurs during one collision between a representative and a physical particle. In these simulations, $\zeta$ is the ratio between physical particles that actually undergo a collision to the total number of physical particles belonging to the swarm of representative particle $i$. The effect of this is that both the dissipated kinetic energy and the time step decreases,

[TABLE]

Another effect of smoothing is that we need to be careful when collisions result in a size change of the representative particle $i$ (coagulation, mass transfer or fragmentation). Since not all particles in swarm $i$ undergo a collision a lot of them have an unchanged size. We resolve this by making use of the Monte Carlo nature of the simulations and draw a random number uniformly distributed between 0 and 1. The probability of selecting one of the particles that has a new radius is $\zeta$. If the random number is >$\zeta$ the collision is treated as a bouncing collision with only energy dissipation as the result. Simulations of collapsing pebble clouds with different values of $\zeta$ show that the formation time of a planetesimal core varies but the overall results are unchanged (Fig. 3, bottom plot). The increased resolution with smaller $\zeta$ causes the energy dissipation to be more evenly (and physically correct) spread out between sub-clouds and it takes longer for the first sub-cloud to contract into a solid core.

2.3 Cloud energies and contraction

In the radially resolved model each sub-cloud corresponds to a representative particle with its own properties (e.g. cloud radius, particle size and speed, etc.). Sub-clouds do not necessarily contract with the same rate (to get into virial equilibrium) and generally do not have the same equilibrium radii, $R_{i,\text{{eq}}}$, as this depends on their individual kinetic and potential energy ($T_{i}$ and $U_{i}$). In virial equilibrium we have

[TABLE]

where $E_{i,\text{{tot}}}$ is the current total energy of sub-cloud $i$. In the old model (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17) the potential was a function of the cloud radius only and the equilibrium radius could be analytically calculated. In the extended model the potential energy of a sub-cloud is more complex; it is the sum of the contributions from all other sub-clouds, $U_{ik}$, and the core (Eqs. 15–24). The expression for $U_{ik}$ is different depending on which sub-cloud, $i$ or $k$, is larger. Now we cannot find $R_{i,\text{{eq}}}$ analytically unless it is smaller than the innermost cloud ($R_{i,\text{{eq}}}<R_{k}$ $\forall k$, Eq. 27) or larger than the outermost cloud ($R_{i,\text{{eq}}}>R_{k}$ $\forall k$, Eq. 29). Otherwise we do not know which sub-clouds $k$ are larger and which are smaller than $i$ at $R_{i,\text{{eq}}}$. This issue is solved by keeping track of the radii and potential energies of all sub-clouds in a sorted array. To, continuously, find the equilibrium radii for sub-clouds we use their total energy and Eq. (6) to get $U_{i,\text{{eq}}}$. Next we compare the equilibrium potential to the sorted radius and potential arrays and find at what cloud size, $R_{i,\text{{eq}}}$, the potential is $U_{i,\text{{eq}}}$. This is done for all sub-clouds every time step as a change in radius of one sub-cloud affects the potentials of all other sub-clouds.

2.4 Collapse limits

Self-gravitating pebble clouds collapse on a short time-scale. The collapse time decreases with increasing planetesimal size but, for cm-sized pebbles, it takes only a couple of orbital periods to form planetesimals of a few km in size (WJJ14). At some cloud mass ($\sim$40 km solid planetesimal radius for cm-sized pebbles) the energy dissipation from inelastic collisions starts to become more rapid than the release of kinetic energy (virialization, Section 2). This leads to a collapse at the free-fall speed. To resolve this in the simulations we limit the contraction speed to the free-fall speed

[TABLE]

where $M_{i,\text{{enc}}}$ is the enclosed mass and $R_{i,0}$ is the radius of latest virial equilibrium of sub-cloud $i$ respectively. The free-fall limit results in, for massive clouds, particles achieving sub-virial velocities and the collapse becomes ‘cold’. One result of this is that even for the most massive pebble clouds collision speeds, in a large part of the collapse, are small enough that collisions do not result in fragmentation and some pebbles survive the entire collapse (WJJ14).

2.5 Effect of gas friction

In the homogeneous cloud model the effect of gas on the collapse was not taken into account. We discussed this approximation in WJ+17 and came to the conclusion that in early parts of the collapse gas drag should matter. In the extended model we include gas drag and limit particle speeds to their terminal speeds, $v_{\text{{t}}}$. The terminal speed for a particle in sub-cloud $i$ can be written as (e.g. Whipple, 1972; Weidenschilling, 1977; Birnstiel et al., 2010)

[TABLE]

where $\tau_{i,\text{{f}}}$ is the friction time of the particle, $F_{i,\text{{g}}}$ is the external force on the particle (gravity), $\rho_{\text{{s}}}$ is the solid density of the particle, $\Sigma_{\text{{d}}}$ and $\Omega_{\text{{K}}}$ are the surface density and Keplerian frequency of the protoplanetary disc at the cloud’s distance to the star. We get Eq. (8) for the terminal speed by assuming that the particles are small (cm-sizes) and experience Epstein drag from the gas. The change in the random speed of a particle in sub-cloud $i$ due to gas drag, $\dot{v}_{i,\text{{d}}}$, is calculated as

[TABLE]

We use the minimum mass solar nebula (Hayashi, 1981) to model the protoplanetary disc containing the collapsing pebble cloud and get $\Sigma_{\text{{d}}}$ to limit both the random speed of particles and the contraction speed of the cloud. Including gas into the simulations can both speed up (particles lose energy faster due to gas drag) and slow down (if the contraction speed is faster than $v_{i,\text{{t}}}$) a collapse.

2.6 Algorithm summary

The algorithm can be summarized in a series of steps done every time step:

Find the time step $\delta t$ and the particles included in the collision (representative particle $i$ and physical particle $k$) from the collision rate matrix $r_{ik}$ (Eqs. 2–3). 2. 2.

Calculate the outcome of the $ik$-collision (Sec. 2.2). Swarm $i$ loses energy $\delta E_{ik}$. The energy lost and time step is decreased from the effect of smoothing (Eqs. 4–5). 3. 3.

Contract all sub-clouds toward virial equilibrium using $\delta t$, potential energy and kinetic energy from the previous time step. To find the equilibrium radius of each sub-cloud we use the sorted potential array (Sec. 2.3). Particle speeds increase due to virialization and decrease due to gas drag (Eq. 9). These contractions are limited by free-fall and terminal speed (Eqs. 7–8). 4. 4.

Increase the mass of the planetesimal core (if any) and decrease the mass of all sub-clouds correspondingly (Eq. 1). 5. 5.

Calculate the new enclosed masses $M_{i,\text{{enc}}}$ (Eqs. 11–14), collision rates $r_{ik}$ and potential energies $U_{i}$ (Eq. 15). 6. 6.

Keep track of free-fall, $v_{i,\text{{ff}}}$, and terminal, $v_{i,\text{{t}}}$, speeds for each sub-cloud so that particles do not move too fast or that the sub-clouds do not collapse too fast. Reset $v_{i,\text{{ff}}}$ if a sub-cloud reaches virial equilibrium.

3 Results

With our simulations we investigate the implementation of the radially resolved collapse model described in Section 2 and how the results differ from the homogeneous cloud model in WJJ14 and WJ+17. We are mainly interested in how the particle size distribution varies with depth inside planetesimals. The shape of the size distribution will affect the packing capabilities and hence the overall density and porosity of the resulting planetesimal. Pebble clouds are formed in gaseous discs around newborn stars and should be embedded in gas. Therefore we also explore the effect of gas on the collapse.

3.1 Initial conditions and simulation setup

For small clouds ($R_{\text{{solid}}}\lesssim$30 km, WJ+17, ) collisions result in bouncing only so pebble sizes remain unchanged. Since the collapse time is a function of pebble size (e.g. WJJ14, ) we are interested in an initial particle size distribution. Investigations of meteorites (e.g. Eisenhour, 1996) suggest that the sizes of chondrules are distributed according to a log-normal or Weibull distribution. We use this as inspiration of the size distribution of pebbles concentrated by the streaming instability and forming self-gravitating clouds. In our simulations we then start with a log-normal particle size distribution (skewed slightly more toward smaller particles),

[TABLE]

where $\mu_{a}$ and $\sigma_{a}$ is the mean and standard deviation of the logarithm of the particle radius $a$. Since we are investigating pebble clouds we set $e^{\mu_{a}}=1$ cm and $\sigma_{a}=0.5$ to study the collapse of a cloud of broadly mm-dm-sized particles.

Initially, the cloud consists a number of homogeneous, spherical, non-rotating, equal-sized sub-clouds of pebbles. It is put at 39 astronomical units from the Sun (the semi-major axis of Pluto, to investigate comet formation) with sizes of each sub-cloud set to the Hill radius of the total mass of all clouds. This choice causes both the density and free-fall time of a pebble cloud to be independent of mass. The radius of the pebbles in each individual sub-cloud is picked randomly from the log-normal size distribution described above.

3.2 Collapsing pebble clouds

The simulations investigate the formation and interiors of comet-sized (5 km) planetesimals from the collapse of gravitationally bound pebble clouds.

Because of higher collision rates (Eqs. 2–3) and stronger coupling to the gas (Eq. 8) smaller particles dissipate their kinetic energy more efficiently than larger pebbles. This causes sub-clouds with different particle sizes to contract at different rates and the density inside the cloud to decrease with distance to the centre (Fig. 2). During the collapse the density increases in the inner parts of the cloud and, at some point, one sub-cloud is the first to reach solid density and will form a planetesimal core. As time goes on, this core accretes particles from the other sub-clouds and grows in mass (Fig. 3, top plot). The increase in core mass decreases the number of remaining pebbles and the density of the outer parts of the cloud which, in turn, results in lower collision rates. In simulations including gas, the particles are subject to gas drag and lose kinetic energy very efficiently even when the number of pebbles is low at the end of the collapse. This keeps densities high (Fig. 2, right plot) and the collapse time short (Fig. 3, top plot). The solid centre of a planetesimal of $R_{\text{{solid}}}=5$ km is formed already after $\sim$25 yr (similar to the free-fall time) and the rest of the material has all been accreted after $\sim$50 yr. For simulations without gas, on the other hand, the particles only lose energy in dissipative collisions resulting in a slower collapse.

An important aim of the collapse simulations is to investigate the porosity inside a pebble-pile planetesimal from the particle packing capabilities. This is done by keeping track of particle size distributions at different depths as material is accreted onto the growing planetesimal core. Fig. 4 shows the particle size distributions at different depths inside resulting planetesimals with radii $R_{\text{{solid}}}=5$ km, one plot for a simulation neglecting gas (left plot) and one with the effect of gas included (right plot). By using these size distributions and a two-dimensional circular version of the drop-and-roll algorithm of Hitti & Bernacki (2013) we produce pebble-pile plots (Fig. 5) to investigate the packing capabilities inside the planetesimals. These plots show the size distribution of pebbles at proper distances from the centre but the porosity in two dimensions does not mirror the porosity in three dimensions. In three dimensions the packing mechanism results in a different porosity. The plots are produced by continuously generating angles randomly between 0 and $2\pi$, drop pebbles toward the centre of the planetesimal from those angles and let them roll to lowest potential (gravity pointing towards the centre). The sizes of the pebbles depends on the current size of the planetesimal and are generated from the size distributions in Fig. 4. Fig. 5 shows that including gas into the simulations has both positive and negative effects on the collapse for pebbles of different sizes. In the case of simulations ignoring gas sub-clouds with small particles contract faster and are the first to form the solid core of the planetesimal. As the collapse goes on, larger and larger pebbles lose their energy to land on the core and the planetesimal grows to have a somewhat layered interior (Fig. 5, left plot). This can also be seen in Fig. 6 (left plot), which shows the particle size vs. cloud size for the sub-clouds at different points in the collapse. The overall trend is that the sub-clouds with large particles take longer time to contract and eventually reach solid density.

In the case of simulations including gas, however, the resulting planetesimal has a different internal structure (Fig. 5, right plot). As can be seen in Fig. 6 (right plot), the sub-clouds with mid-sized ($\sim$1 cm) pebbles collapse the fastest and form the core. Both sub-clouds with smaller and larger particles, on the other hand, take longer to reach solid density and, mixed together, form the outer layers. This is caused by the fact that gas not only dissipate energy from the pebbles but also limits the collapse speed to the terminal speed (Eq. 8) of the particles. The sub-clouds with the smallest particles lose their kinetic energy quickly by gas drag but cannot contract faster than the terminal speed. The larger particles have a higher terminal speed ($v_{\text{{i,t}}}\propto a_{\text{{i}}}$) and the contraction is less limited by the gas. However, the gas drag on the larger particles is smaller and they lose energy less efficiently which eventually result in a slower collapse. The optimal particle size is somewhere in between ($\sim$1 cm) where the terminal speed is high enough to not hinder the collapse too much but still small enough that the gas drag is efficient.

One difference between the radially resolved (this paper) and the old homogeneous model (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17) is that fragmenting pebble-pebble collisions also occur in the formation of low-mass planetesimals. The surface layer of a resulting planetesimal contains small particles (pebble fragments) which can be seen in the particle size distributions of the outermost layer in Fig. 4. As discussed in Section 2.4 the collapse becomes ‘cold’ (sub-virial particle speeds) if the energy dissipation is too rapid. This helps pebbles to survive collisions in the end of the collapse where the virial speed needs to be high to support the increasing gravity (Eq. 6). The last phase of a cloud collapse consists of a solid planetesimal core surrounded by a small fraction of the total mass as a cloud with low density. The low density results in slower energy dissipation and the possibility to get closer to virial equilibrium after each collision. An increase in terminal speed (Eq. 8) as cloud contracts allows for higher particle speeds even when the effect of gas is included. In turn, the collision speeds reach values high enough for fragmentation and small pebble fragments will land on the surface of the planetesimal. This does not happen in the old homogeneous model (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17) because there the sub-clouds collapse together resulting in efficient energy dissipation in the end of the collapse (a core with surrounding low-density cloud never forms).

3.3 Pebble-pile comets

The formation of comets (and other km-sized planetesimals) through the collapse of self-gravitating pebble clouds results, as in (\al@wahlberg14,wahlberg17; \al@wahlberg14,wahlberg17), in porous pebble-piles. In low-mass pebble clouds collision speeds only reach values high enough for silicate dust aggregates to fragment in the very end of the collapse (Fig. 4). Ice aggregate pebbles would not fragment at any point in the collapse to small planetesimals (WJ+17). Our simulations show that the environment of the initial pebble cloud is important to the internal structure of the resulting planetesimal. If the disc gas surrounding the pebbles is taken into account, the collapse process is different. Gas increases energy dissipation through drag but also limits the contraction speed of pebbles (Sec. 2.4). This causes a comet formed in the outer regions of the Solar System to be a pebble-pile with mid-sized ($\sim$1 cm) pebbles in the centre of the comet and surface layers made of a mixture of larger and smaller pebbles (Fig. 5, right plot). A mixture of particle sizes would lead to more efficient packing and lower porosity. Without disc gas the comet has an ‘onion’-like interior with shells of pebbles decreasing in size with increasing depth (Fig. 5, left plot). The shell structure observed for each individual lobe of 67P/Churyumov–Gerasimenko (Massironi et al., 2015) suggests the lobes could have formed from the contraction of a pebble cloud and that activity erodes into layers of different pebble size distributions (and hence porosities). Measurements of a permittivity gradient down to $\sim$100 m inside 67P (Ciarletti et al., 2015) implies a variation comet bulk properties. The permittivity decreases with increasing depth suggesting higher porosity and/or decrease in dust-to-ice ratio further inside the comet.

4 Conclusions

In this paper we have expanded the Wahlberg Jansson et al. (2017) model of planetesimal formation through the collapse of gravitationally bound pebble clouds inside protoplanetary discs. The collapse results from energy dissipation in inelastic pebble-pebble collisions inside the clouds. Our previous models have assumed that particles were spread homogeneously inside the cloud. Our new model takes into account the presence of the gaseous protoplanetary disc and the radial structure of the resulting planetesimal. To investigate the effects of disc gas we run two sets of simulations, one with and one without the effects caused by gas (Sec. 2.4).

As suggested in the analysis of Wahlberg Jansson & Johansen (2014) the collapse process is a strong function of particle sizes. The total surface area is larger for small particles than for an equal mass in larger particles and this results in higher collision rates and more efficient energy dissipation. Small particles are also more affected by gas (Eq. 8), increasing the energy dissipation further. This causes the particle size distribution inside a planetesimal to be a function of depth. Most of the largest pebbles are found in the outer layers of the planetesimal (Fig. 4). In simulations ignoring gas this results in a planetesimal with an ‘onion’-like structure. Including gas, however, leads to a planetesimal interior with mid-sized pebbles ($\sim$1 cm) in the centre and smaller and larger particles in the outer layers. This is caused by the fact that gas not only provides drag (stronger for small particles) but also limits the contraction speed to the terminal speed of the particles (Eq. 8). The particle size distribution inside the planetesimal is also affected by the evolution of the density of individual pebbles.

An important effect of resolving the radial structure is that a core forms early in the collapse phase and grows by accreting the surrounding cloud. Without gas, energy is dissipated only by particle collisions and the collapse time-scale is long (Fig. 3, top plot). For the gas model, gas drag leads to continuous energy dissipation even at low densities towards the end of the collapse. At the end of the collapse phase most of the mass is in the core and the remaining cloud has a low density. The low density leads to low collision rates and the possibility to get back into virial equilibrium, leading to higher collision speeds which, in turn, cause fragmenting collisions. This results in a broader particle size distribution in the outer layers of the planetesimal Fig. 4. In previous models the collapse is cold (sub-virial velocities because of free-fall limit) so, for small planetesimals, bouncing is the only outcome of collisions.

Our results suggest that comets formed through the collapse of pebble cloud are pebble-piles with a porous interior. The sizes of pebbles will vary with depth inside the comet (Fig. 4). The resulting comet structure differs if the cloud collapsed embedded in gas or not Fig. 5. Our simulations assume pebbles are silicate dust aggregates which result in possibility for fragmenting collisions at the end of the collapse. If the pebbles are made of ice, however, the collision speeds would never reach values high enough for fragmentation (Wahlberg Jansson et al., 2017) and the range of particle sizes in the surface layer of the comet would be narrow.

Our radially resolved collapse model still lacks physical properties, such as rotation. A non-zero initial angular momentum of a pebble cloud undergoing a collapse could result in binary objects or discs forming around new planetesimals. Through N-body simulations Nesvorný et al. (2010) found formation of binary planetesimals in gravitational instability of pebble concentrations. This planetesimal formation mechanism matches observations of trans-Neptunian objects. Both a high fraction of binary objects and that components in a binary have the same composition have been observed (Noll et al., 2008; Benecchi et al., 2009). The next step in our investigations is to increase the dimensionality further and investigate the effect of rotation. discs should form around planetesimals which we plan to investigate with hydrodynamical simulations and present in a subsequent paper.

Acknowledgements

KWJ and AJ were supported by the European Research Council under ERC Starting Grant agreement 278675-PEBBLE2PLANET. AJ was also supported by the Swedish Research Council (grant 2010-3710) and by the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017 and 2014.0048). KWJ and AJ acknowledges the support from the Royal Physiographic Society in Lund for grants to purchase computer hardware to run the simulations on. We also thank the referee Kees Dullemond for insightful comments that helped improve the manuscript.

Appendix A Enclosed mass and potential energy of sub-cloud $i$ in the radially resolved pebble cloud model

In the radially resolved model described in the main text the enclosed mass for each sub-cloud $i$, $M_{i,\text{{enc}}}$, can be written as

[TABLE]

where $m_{i}$ and $R_{i}$ is the mass and radius of sub-cloud $i$ and $M_{\text{{c}}}$ is the mass of any solid core. We get the total potential energy for sub-cloud $i$, $U_{i}$, from Poisson’s equation,

[TABLE]

where $\rho_{i}$ is the mass density inside sub-cloud $i$. Compared to the homogeneous model in WJJ14 and WJ+17 the equilibrium radii in the extended model cannot be found analytically, except for the smallest and largest sub-cloud. As described in the main text (Sec. 2.3) this is solved numerically by keeping track of the potential energies of all sub-clouds in a sorted array. For a cloud in virial equilibrium we have

[TABLE]

where $E_{i,\text{{tot}}}$ is the current total energy of sub-cloud $i$. If $i$ has an equilibrium radius within all other sub-clouds ($R_{i,\text{{eq}}}<R_{k}$, $\forall k$) it can be calculated through Eq. (15) and Eq. (25),

[TABLE]

and solve for $R_{i,\text{{eq}}}$. Similarly we can find $R_{i,\text{{eq}}}$ if $i$ is larger than the radius of all other sub-clouds, ($R_{i,\text{{eq}}}>R_{k}$, $\forall k$),

[TABLE]

and solve for $R_{\text{{i,eq}}}$.