© ESO 2021

1 Introduction

The formation mechanism of planetesimals and comets remains uncertain. While micrometer grains can grow through coagulation to roughly centimeter sized pebbles (Dominik & Tielens 1997; Birnstiel et al. 2012), growth stagnates due to the bouncing barrier (Zsom et al. 2010) and the fragmentation barrier (Blum & Münch 1993). However, if growth to meter sized boulders is obtained by coagulation, these objects drift toward the star within several hundred orbital timescales due to the nebular gas (Weidenschilling 1977a; Nakagawa et al. 1986). It was long thought that the increased sticking properties of water ice could circumvent the bouncing barrier (Wada et al. 2009; Gundlach et al. 2011), though recent lab experiments have shown that this advantage only holds in a rather narrow disk temperature range (Musiolik & Wurm 2019) and that relative velocities quickly result in fragmentation (Blum & Wurm 2008).

The currently favored growth model is the gravitational collapse of a pebble cloud due to highly concentrated clumps of solids (Johansen et al. 2007, 2009) induced by the streaming instability (SI; Youdin & Goodman 2005). The high solid-to-gas ratios needed to trigger SI are not straightforward to achieve in typical protoplanetary disks. Studies do however show that so-called pressure bumps lead to over-densities of solids (Whipple 1972; Brauer et al. 2008; Drążkowska et al. 2013). Furthermore the re-condensation of icy pebbles just outside the snowline can lead to a pileup of solids, enhancing pebble surface densities by a factor of five or more (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017).

A limited spatial resolution in large-scale SI simulations creates difficulties in following the final collapse phase of these formed clumps. This phase is vital in understanding the final structure of the solid core and how this compares to observations of planetary bodies. N-body simulations of a collapsing spherical cloud of pebbles with significant rotation produce similar mass binaries, such as Pluto and Charon (Nesvorny et al. 2010a). Also, recent large-scale hydrodynamical simulations match the 80 percent prograde binary rotation direction observed in KBOs (Nesvorný et al. 2019).

Detailed observational data from comets obtained by the ESA Rosetta probe provide comparison material for the pebble cloud collapse model as a potential comet formation mechanism. Comets have a porosity of up to 80 percent, are several kilometers in size, have a low bulk density of roughly 0.5 g cm⁻³ (Groussin et al. 2019), and appear to have a layered pebble-pile structure throughout their nucleus ranging from one to several tens of millimeter (Poulet et al. 2016). Several numerical studies have been performed to explain these properties by postulating the formation of comets through the collapse of a pebble cloud (Blum et al. 2017).

Pebbles dissipate kinetic energy through collisions and gas drag, which leads to cloud contraction, as shown in pioneering work by Wahlberg Jansson & Johansen (2014, 2017). Depending on the total cloud mass, the initial size distribution can be altered significantly during collapse due to collisions that alter pebble mass, as well as varying terminal velocities for different pebble Stokes numbers. Even bouncing collisions can affect the initial imposed size distribution through pebble compaction,which aids the SI if the compaction timescale is faster than the radial drift timescale toward the star (Lorek et al. 2016, 2018).

We takea more direct approach by solving the dynamical equation of motion of individual pebbles (swarms) that are subject to gas drag, mutual collisions, and gravity. In particular, we are interested in the collapse timescale and the final radial pebble size distribution over the collapsed core. We compare our findings with earlier work and observational data to increase understanding of the formation and final structure of comets and planetesimals.

The paper is structured as follows. In Sect. 2, we present the general setup of our cloud collapse model. In Sect. 3, we present the results of our simulations followed by the discussion in Sect. 4 and a summary and conclusions in Sect. 5.

2 Model and setup

2.1 Dynamical cloud evolution

To perform the numerical models, we developed the code implode¹. We follow the collapse of a cloud of self-gravitating and colliding pebbles resulting from the SI (Youdin & Goodman 2005). We consider a spherically symmetric cloud of pebbles² with a distribution of masses, f(m). The total cloud mass in pebbles M_t corresponds to the mass of a solid planetesimal core of radius R_c and density ρ_• = 1 g cm⁻³. The pebbles are distributed evenly over the cloud’s Hill sphere, which is defined as: $$R_{\text{H}}=r_0{\left(\frac{M_{\text{t}}}{3M_{\star }}\right)}^{\frac{1}{3}},$$(1)

where r₀ is the orbital distance from the star with mass M_⋆. The pebbles interact through mutual gravity, gas drag and mutual collisions. We extended our model to 3D to incorporate the effect of initial dispersion in the angular directions of the spherical cloud. To focus on the analysis of the influence of collisions and gas drag, we neglect shearing effects in our model. This will be incorporated in future studies to investigate the effect of increasing complexity.

The spherical symmetry imposed on our cloud allows us to model the gravitational interactions with the shell theorem: Pebbles only feel the gravitational pull of the mass distribution situated below them as if this mass were concentrated at the spherical center as a point source. Thus the gravitational acceleration of a pebble at radius r_p is calculated from the mass situated at r < r_p with respect to the cloud center of mass: $$f_{\text{g}}=-\frac{G{M}_{\text{enc}}}{r_p^3}{r}_{\text{p}},$$(2)

where M_enc denotes all mass that is closer to the center of mass than the corresponding pebble at position vector r_p in the cloud. We took the effect of gas into account with a simple drag prescription of the form: $$f_{\text{d}}=-\frac{1}{t_{\text{s}}}{v}_{\text{p}},$$(3)

where v_p is the velocity vector of the pebbles. For pebbles of radius s and internal density ρ_• the stopping time is given by (Whipple 1972): $$t_{\text{s}}=\{\begin{array}{cc}\frac{{\rho }_{•}s}{{\rho }_{\text{g}}{v}_{\text{th}}}& \text{Epsteinregime}:s<\frac{9}{4}{l}_{\text{mfp}}\\ \frac{2{\rho }_{•}{s}^2}{9{\eta }_d}& \text{Stokesregime}:s\ge \frac{9}{4}{l}_{\text{mfp}}\end{array},$$(4)

where v_th is the thermal speed, ρ_g the mid-plane gas density, l_mfp the molecular mean free path, and η_d the kinematic viscosity of the gas. The Epstein regime is the relevant regime for our range of pebble sizes since for 10 au the molecular mean free path is several meters. The Stokes number (St) relates the stopping time to one orbital timescale: $$\text{St}=t_{\text{s}}{\Omega }_0,$$(5)

with Ω₀ the Kepler orbital frequency. The equation of motion governing the dynamics of these pebbles is given by: $$\frac{\text{d}{v}_{\text{p}}}{\text{d}t}=-\frac{G{M}_{\text{enc}}}{r_{\text{p}}^3}{r}_{\text{p}}-\frac{v_{\text{p}}}{t_{\text{s}}}.$$(6)

Fig. 1 2D sketch of the spherical cloud collapse model with the center of mass (CM) in theorigin. Pebble swarms (blue circles) are initiated over a Hill sphere with spherical velocities and positions. During evolution pebble swarms feel both gravity fg and gas drag fd. The dashed shells indicate the zones for collision evolution. As the pebbles settle, zones are being rebuilt such that the number of swarms per zone never falls below a desired value Nz.

2.2 The representative particle approach

The number of physical pebbles in the collapsing cloud is on the order of 10¹⁷ in our models. To be able to model them numerically, we model pebble collisions and advection using the representative particle approach (Gillespie 1975; Zsom & Dullemond 2008) with much lower number of N pebble swarms³, where N ≪ N_phys. Every pebble swarm represents N_i physical pebbles with equal properties. The mass M_swarm = M_t∕N of each swarm is equal and remains constant at all times.

We use an adaptive grid developed by Drążkowska et al. (2013) to resolve the collisions locally in the cloud. We assume the cloud to be spherically symmetric. The natural way to subdivide local collision zones is through spherical shells over the cloud domain. The number of pebble swarms per spherical shell is kept above a desired value N_z by rebuilding the shells if needed, to ensure enough resolution to perform collisions in every shell. In Fig. 1, we provide a sketch of the cloud collapse model. Collisions are then modeled with a Monte Carlo algorithm based on the total “particle in a box” collision rate in a zone $$R={\sum }_i^{}{\sum }_k^{}{R}_{ik},$$(7)

where R_ik = n_kσ_ikΔv_ik is the collision rate between a representative pebble from swarm i and a nonrepresentative pebble of swarm k.

The collision time step is then chosen from the collision rate by drawing a random number U ∈ [0, 1) as: $$\delta t_c=-\frac{1}{R}\ln (U).$$(8)

After performing the collisions over the time step, the matrix of collision rates is updated accordingly. Since it is unlikely that representative pebble i collides with any pebble from swarm k we follow, we only update the rate of representative pebble i. In the case of mass change of pebble i during a collision the physical pebble number represented by swarm i changes to N_i = M_swarm∕m_p where m_p is the pebble massand lower in case of fragmentation, and higher in case of coagulation.

For a fragmentation threshold velocity below 10 m s⁻¹, or if bouncing collisions stop growth a very low Stokes numbers, the SI will not be triggered (Drążkowska & Dullemond 2014; Drążkowska et al. 2016). To meet the context of SI simulations, we limit ourselves in the possible collisional outcome, as a function of collisions velocity to sticking: Δv ≤ 10 m s⁻¹ and fragmenting: Δv > 10 m s⁻¹ (Brauer et al. 2008). After a sticking collision the mass m_i of the representative particle is updated to ${m^{\prime }}_i=m_i+m_k$ and we use conservation of momentum to calculate its resulting velocity. In case of fragmentation, we assume thatthe mass of the original representative particle is distributed according to the power-law n(m)∝ m^−11∕6, consistent with the MRN size distribution (Mathis et al. 1977). We choose the new mass of the representative particle randomly from these fragments.

Fig. 2 Overview of the algorithm starting from cloud initialization at t = 0.

2.3 Initial conditions

We initiate a cloud of pebbles at a orbital distances of 10 and 39 au around a solar mass star, respectively. The pebbles are given positions and velocities (3D) (r_p, ϕ_p, θ_p, v_r, v_ϕ, v_θ) where we have adopted a spherical coordinate system. In the radial direction r_p, pebbles are placed such that the volume density is constant over the cloud. In the angular directions ϕ_p, θ_p, the positions are uniformly randomized over the angle domains ϕ ∈ [0, 2π] and θ ∈ [arccos(−1), arccos(1)], respectively.

The SI in general relies on pebbles between minimum and maximum Stokes number St_min ~ 2 × 10⁻³, St_max ~ 10⁻¹ (Bai & Stone 2010; Yang et al. 2017). We use an initial MRN size distribution between these Stokes numbers (Mathis et al. 1977). For the chosen size distribution most of the cloud mass is found in the largest pebble sizes, while the lower sizes contain the most pebbles in number. Pebbles in these ranges of Stokes numbers are influenced by the gas on a timescale that is short compared to the cloud collapse timescale. We therefore initiate pebbles with their corresponding terminal velocity in the radial direction (see Appendix A for details). For the angular directions (v_θ, v_ϕ), following Wahlberg Jansson & Johansen (2014), we initiate pebbles with a Maxwellian velocity dispersion: $$\text{d}P(\Delta v)=\frac{1}{2\sqrt{\pi }}\frac{\Delta v^2}{{\sigma }^3}{e}^{-\Delta v^2/4{\sigma }^2}\text{d}(\Delta v),$$(9)

with the dispersion $\sigma =\sqrt{2K_0/M_t}$. The random velocities are extracted from the initial potential U₀ and kinetic energy K₀ of the cloud, assuming initially virial equilibrium. A schematic overview of the cloud evolution is given in Fig. 2. The gas in the disk is modeled according to the minimum mass solar nebula (Weidenschilling 1977b; Hayashi et al. 1985) with power law expressions for the gas temperature and surface density, respectively: $$T(r_0)=170\text{K}{\left(\frac{r_0}{1\text{au}}\right)}^{-1/2},$$(10) $$\Sigma (r_0)=1700\text{g}{\text{cm}}^{\text{-2}}{\left(\frac{r_0}{1\text{au}}\right)}^{-3/2}.$$(11)

Simulations of the SI show that the planetesimals formed have a typical radius on the order of 50–100 km (Johansen et al. 2009; Schäfer et al. 2017). Kilometer-sized comets are unlikely to form via the SI in turbulent environments as clumps are easily broken up again. An increase in resolution of the SI simulations shows, however, that the smallest clumps decrease in size (Johansen et al. 2011; Simon et al. 2016). Another reason that km sized object could survive is the fragmentation of a single clump into binary or satellite component (Nesvorný et al. 2010b, 2019). We consider three different cloud masses spanning the lower uncertain constraint on the core radius R_c = 1 km and the more certain upper constraints R_c = 10 km and R_c = 100 km. We locate the clouds at orbital distance 39 au (Kuiper belt distance) and 10 au (Saturn distance). We show the Hill radii compared to gas scaleheight H_g and the core masses M_Ceres in Table 1. An additional table with the meaning of symbols is shown in Table D.1.

2.4 Relative collision speeds

Modeling collisions accurately in this system is nontrivial because the velocities in the radial direction are systematic. As the cloud collapses under the forces of gravity, pebbles are accelerated toward the center, with velocities moderated by the gas drag. Whether pebbles are moving with terminal (due to friction) velocities or in free fall, these velocities strongly depend on location.

In a representative particle approach, we bin the radial structure of the cloud and consider collisions between pebbles in a radial bin. We consider collisions between randomly selected pebbles in this bin, looking at relative velocities that may cause collisions. Due to the systematic nature of the radial motions, and in particular the radial gradient in velocities, the relative velocity of any two pebbles will be overestimated unless the pebbles are located at exactly the same distance from the center of mass. Because of this effect, we can expect that the results would significantly depend on the resolution in the model. To avoid this, we assume that the radial velocity of each pebble, for the properties of the collision, is the one the pebble would have when located exactly in the middle of each cell. In addition to that, we still consider random velocities in the θ and ϕ direction. These velocities are initially set from the initial conditions, are damped by the interaction with the gas, and are fed by the outcome of collisions. In Appendix C we show that this approach leads to robust convergence of the results as a function of resolution.

2.5 Accretion conditions

Another numerical challenge of this system is that the core forming in the center is orders of magnitudes smaller than the starting size of the cloud. It is important to test if pebbles approaching the center of mass will get stuck there due tocollisions right away, or if they will pass through the center and start a damped oscillation before finally settling toward the core. We test this by checking at what point the center of the collapsing cloud becomes so dense that no pebbles will be able to pass through without suffering collisions. Basically, this is the same as asking, when the core will become optically thick. For our computation, we will assume that if the first pebbles arriving in the core become optically thick, a core has formed. After that, any pebbles arriving within a set distance of the center of mass will settle onto the core in the sequence of arrival. The rapid formation of a massive core was also observed by Wahlberg Jansson & Johansen (2017).

In all models, we are tracking the optical depth of one percent of the total cloud mass M_t to ensure that this assumption is reasonable. In Appendix C, we demonstrate that this condition is matched for our fiducial model. If the innermost one percent of optically thick pebble swarms has settled, the cloud radius at this distance is adjusted to be the accretion radius for the outer swarms as r_acc = r_τ>1. Here r_τ>1 is the position where the first one percent inner mass reached τ_in > 1. Pebbles are therefore accreted if they meet the condition: $$r_p-r_{\text{acc}}<0.$$(12)

with which we settle pebbles in the core with uniform density, leading to growth in core radius according to $r_{\text{c}}={(3N_{\text{sett}}{m}_{\text{p}}/4\pi {\rho }_{•})}^{1/3}$, with N_sett the amount of pebbles settled. We also keep track of the sequence with which pebbles settle in time to preserve information about the final size distribution of pebbles over the radial extend of the core layers.

2.6 Numerical scheme

We integrate the equation of motion of the pebble trajectories while we simultaneously resolve the mutual pebble collisions using the representative particle approach. The equation of motion is integrated in 3D with a Runge-Kutta Fehlberg variable step scheme (Fehlberg 1969) with an error tolerance tol = 10⁻⁷. The time step following from the RKF45 solver, δt_EOM is compared with the longest time step we can afford to include collisional evolution reasonably well, δt_c, which is calculated as a minimum of timesteps reported by each radial zone: the average timestep between two consecutive collisions multiplied by the number of particles in the zone. We pick the minimum time step from this to ensure we resolve both the dynamical and collision part of the cloud evolution: $$\Delta t=\min [\delta t_{\text{EOM}},\delta t_{\text{c}}].$$(13)

If pebbles are within the accretion radius r_acc we place them in the core accordingly. The N_sett pebbles are considered as settled do not participate in collisions and in the further cloud evolution anymore.

3 Results

In Figs. 3 and 4, we present results for the cloud located at 39 au (with R_c = 1 km the fiducial value) and 10 au orbital distance from the central star, respectively. An overview of the parameter study is given in Table 2. The total number of representative pebbles is N = 10⁴. The core density and internal pebble density are 1 g cm⁻³ and the initial size distribution for all three cases corresponds to [St_min = 2 × 10⁻³, St_max = 10⁻¹].

Starting with the 39 au case, we show the average number of collisions per representative pebble $\overline{N}=N_{\text{coll}}/N$, with N_coll in Fig. 3, top left panel. For increasing core radius R_c, $\overline{N}$ increases due to the higher number density in smaller pebbles. An estimate of the increase rate of $\overline{N}$ for larger core radius is given in Appendix B. Starting from t = 0, the average collisions number $\overline{N}$ increases and approaches a constant value after some time. This is explained as follows. The most extreme difference in radial velocities is found between the largest and the smallest pebbles. The local collision rates in the radial direction are fed predominantly by this extreme case. We typically observe that the most massive pebbles form the first core, leading to a large drop in Δv_r for the smaller active pebbles. By the time when $\overline{N}$ is constant, the core has formed from the inner one percent of cloud mass with τ_in ~ 1. This sets the accretion radius for the active remaining pebbles in the cloud. Collisions that would normally occur in the over-dense region within the accretion radius are now excluded. While this could lead to some collisional processing in the collapsing shells, the high optical depth guarantees that pebbles cannot change sequence within r_acc.

The top right panel of Fig. 3 shows E = K∕|U|, the ratio of the total kinetic energy to the total absolute value of the potential energy of the pebbles in the cloud. We observe that the system never reaches anything close to the virial equilibrium, which would correspond to E = 0.5, except for the rapid interval where freefall restricts addition of further energy at the peak of E (top right panel of Fig. 3, global peak). Instead, the formation of a massive core happens before pebbles can cross the central region of the cloud.This prevents pebbles from oscillating through the center of the cloud and causes the kinetic energy of the system to drop to zeroif all pebbles have settled.

The evolution of E is explained as follows. Starting at t = 0, the kinetic energy of the individual pebbles is fed by the virial assumption for initial dispersion given in Eq. (9). As time evolves, the most massive pebbles typically form the first core characterized by the sharp global maximum in E. The fall times of the massive pebbles are approximately equal in the narrow bigger size range (Fig. A.1) and they are only weakly coupled to the gas. The sharp global maximum in E is a combination of the potential energy efficiently being converted to kinetic energy due to their high terminal velocities and a collective, near instantaneous, arrival at the center of mass. Shortly after the maximum, these massive pebbles have formed the core and no longer contribute to the E evolution, explaining the steep drop in E. The difference in fall times for smaller pebbles increases rapidly with decreasing size since $\frac{\text{d}{t}_f}{\text{dSt}}\propto -1/{\text{St}}^2$. This causes a much more gradual decrease in E. If pebbles have reached the accretion radius, they are frozen in the core. The total energy therefore declines with respect to the energy at t = 0 since we exclude them from the further energy evolution of active pebbles. This explains the fluctuation in total energy during the decline to E = 0, every time apebble is settled. In Appendix C, we present a resolution study for different representative pebbles. The most telling parameter is the average number of collisions experienced by a pebble, and we can see that, for sufficient resolution, $\overline{N}$ remains the same for increasing N.

The core mass fraction (CMF; Fig. 3, bottom left panel) indicates the fraction of mass in pebbles that has settled with respect to the total cloud mass (bottom-left panel). The steep vertical increase in the CMF can be understood in the same way as the global peak in E. The massive pebbles all fall to form the core together, leading to a sudden increase in the CMF. The more gradual infall of the smaller pebbles leads to a smoother and more gradual increase in the last few percent of the CMF referring to the same explanation as for the behaviour in E after the global peak. For massive clouds that are optically thick in the inner region already at t = 0, the CMF may be nonzero already from the start since we form the core at t = 0 in this case. This is only observed for the 10 au case for R_c = 100 km (Fig. 4, bottom left panel). The initial and final size distribution is shown in Fig. 3, the bottom right panel. For low mass clouds (R_c = 1 km) the average collisions per pebble $\overline{N}~10$ shows that pebbles are hardly affected by collisions, even more so since collisions are dominated by larger pebbles colliding with smaller ones. This reflects back in the final size distribution being almost identical to the initial one. As cloud mass increases, collisions become more important and the final size distribution is shifted increasingly toward higher Stokes numbers. This shows that the massive pebbles grow larger by sweeping up the smaller ones through coagulation since the small pebbles from the initial distribution are depleted in the final one. It is worth noting that we are not in the classical situation of coagulation happening over orders of magnitude in size with a narrow size distribution, in which even a low number of collisions can cause significant growth. Rather, we have large pebbles collecting some smaller ones. So even in the case of collecting pebbles of the same size as the original one, $\overline{N}$ collisions only increase the mass by at most a factor of $\overline{N}$. $\overline{N}~10$ means a maximum change of size, and the proportional Stokes number, by a factor of 10^1∕3 ≈ 2. This is consistent with the Stokes number distribution presented in the bottom right panel of Fig. 3.

For the computations at 10 au, we show the same parameters as the 39 au case in Fig. 4. The most noteworthy change in results is a shorter fall timescale for the cloud due to decrease in cloud Hill radius by a factor of four. The same mass is now initially distributed over a smaller region. Average collisions per pebble inevitably increase with respect to the 39 au case (Fig. 4 top left panel). This leads to a more rapid increase of the optical depth in the central region. For the 100 km case, the inner active core mass already becomes optically thick at t = 0. The CMF increase is therefore, after an initial jump, more gradual already at t = 0 since the full mixture of the initial size distribution contributes immediately by falling through the accretion radius r_acc (Fig. 4, bottom left panel). An additional overview of the main parameters used in the simulations is given in Table 2.

Fig. 3 Simulation results for three different cloud masses corresponding to core radius Rc = [1, 10, 100] km at 39 au from the central star. Top left panel: average collisions per pebble $\overline{N}$. Top right panel: ratio of kinetic to potential energy E = K/|U|, bottom left panel: CMF in time. Bottom right panel: mass weighted initial and final pebble Stokes number distribution. The initial St distribution (red curve) is the same for all three cloud masses.

3.1 Fragmentation events

We rarely observe fragmentation happening in the simulations, even for very massive clouds in which pebbles should, in principle, reach terminal velocities exceeding the fragmentation threshold Δv = 10 m s⁻¹. To find out why this is the case, we look at the maximum terminal velocities that pebbles can reach before arriving at the accretion radius. We estimate the accretion radius where the optical depth of the first inner one percent of the cloud mass M_t exceeds unity (see Appendix C for details). $${\tau }_{\text{in}}=\frac{N\sigma }{4\pi r^2},$$(14)

with N = 0.01M_t∕m_p the physical number density of pebbles contained in the inner one percent of total cloud mass, σ = πs² their physical cross-section and 4πr² the surface of the sphere enclosing the 0.01M_t in pebbles. If τ_in > 1, the contracting core is optically thick and becomes impenetrable for outer mass. Solving for the accretion radius r_acc in Eq. (14) gives $$r_{\text{acc}}={\left(\frac{0.01R_{\text{c}}}{4s}\right)}^{1/2}.$$(15)

Filling in this expression for the terminal velocity of a pebble at this particular radius $v_{\text{t}}=G0.01M_t\text{St}{\Omega }_0^{-1}/r_{\text{acc}}^2$ yields a critical terminal velocity v_t,* at the accretion radius of: $$v_{t,*}~60\text{m}{\text{s}}^{\text{-1}}{\text{St}}^2.$$(16)

This expression is independent of disk radius for the specific power law n = −3∕2 of the gas surface density (Eq. (11)) and only depends on the square of the Stokes number. For other power law profiles of the gas surface density there will only be a weak dependence.

To summarize this result: fragmenting collisions between pebbles are difficult to achieve due to optical depth of the cloud exceeding unity before the collapse has been finalized (Appendix C, Fig. A.3, bottom right panel). The cloud becomes optically thick at cloud radii where the terminal velocity needed for fragmentation cannot be reached. We expect (even though we are unable to track these collisions for numerical reasons) that pebble collision rates beyond this region prevent pebbles from reaching high speeds due to rapid damping. We propose from these findings that fragmentation events in collapsing pebble clouds are rare and comets are composed from primordial building blocks. Indeed this result is consistent with models of the surface and interior of comet 67P that reveal an active surface layer of ~1 m and primordial pebbles below that (Capria et al. 2017). This strengthens the claim that small solar system objects are in general primordial in nature (Bottke et al. 2005).

Fig. 4 Same results as shown in Fig. 3 but now for an orbital distance of 10 au for N = 104 in all three cases. The most notable changes are a significantly shorter collapse time, increasing number of collisions due to smaller Hill radii, and a gradually increasing CMF. Growth can be significant as a result of sticking collisions for 100 km core radius.

3.2 Finalcore structure

In Fig. 5, we show 2D XY slices through the center of the final core for the fiducial model with R_c = 1 km, R_c = 10 km, R_c = 100 km and an extreme case of R_c = 500 km at 39 au. In Fig. 6, we present the contribution of the radial core layers to the total core mass for 1 and 100 km core radii, at 39 and 10 au. In general, the resulting core structure for low mass clouds (R_c = 1, 10, 100 km for 39 au and R_c = 1, 10 km for 10 au) can be explained as follows. The sequence of arrival of pebbles is dominated by their corresponding terminal velocities. Pebbles are found further from the center of mass of the core for decreasing terminal velocity. The resulting core structure is characterized by massive pebbles in the center, followed by increasingly smaller pebbles toward the core surface (Fig. 5, top left, top and bottom right panel).

The massive pebble swarms form the core, and while they dominate the mass, they are fewer in number. This explains that the massive inner core region is only a small mass contribution with respect to the total core mass (Fig. 6, the top and bottom left panels). If we look closer toward the core surface, we observe that contributions to the total mass become higher. The main reason for this is that the pebble number density increases for decreasing Stokes number. It then follows that more and more pebbles fall at approximately the same time to the core. The total amount of pebbles in one pixel of width Δr and height Δs increases toward the surface since the mass of one swarm is constant.

In general, in all cases in Fig. 6 there is a region with pebbles larger than the pebbles in the core while they have ended up further outwards. These are pebbles that grew through coagulation while the largest pebbles from the initial size distribution already formed the core. The sequence of increasingly larger pebbles with increasing core depth is not fully observed in the 100 km (10 au) and 500 km (39 au) core (Fig. 5 bottom panels, Fig. 6 bottom right panel). For the largerst core mass M_t at a given orbital distance, the inner region becomes optically thick at an earlier stage due to τ_in ∝ M_t. This sets the accretion radius at an increasingly higher fraction of the Hill radius for larger M_t (r_acc ~ 0.2R_H for R_c = 500 km).

The observed accretion sequence of pebbles in the massive cores reflects the initial placement of the pebbles. An extreme example is the abovementioned run with 500 km planetesimal at 39 au, for which the inner part of the pebble cloud is optically thick already at t = 0. All the pebbles within the accretion radius r_acc ~ 0.2R_H are directly settled including the small pebbles. Additionally, since pebbles now accrete at r_acc ~ 0.2R_H, there is less time to differentiate the fall timescales of the small and large pebbles resulting in a more randomized accretion sequence for the outer core too. This observation implies that a significant core might already form during the SI phase for massive clouds, before the gravitational collapse of the full cloud is even starting.

An important aspect to mention is our choice of initial cloud radius. As mentioned by Wahlberg Jansson & Johansen (2017), the choice of distributing the cloud mass over its corresponding Hill radius leads to self-similar collapse for given orbital distance r₀. The reason for this is that free-fall timescales as well as terminal-velocity fall timescales are independent of cloud density and mass. We therefore predict that comets and planetesimals are composed of a solid core formed by primordial massive pebbles. The outer layers are constructed of increasingly smaller pebbles toward the surface of the body. Especially for comets the outer layers would be stripped down easily over tens of eccentric orbits around the central star, revealing the inner core of primordial centimeter sized pebbles (Pajola et al. 2017; Arakawa & Ohno 2020).

Fig. 5 2D slices in the XY plane of the final core formed from the collapse on 39 au normalized in units of Rc. The z-dimension has been flattened to one percent of the core radius. Pebbles are indicated with the circles and scale from smallest circle to largest circlewith Stmin, Stmax resp. Top left panel: core structure for Rc = 1 km. The inner core is composed of primarily Stmax pebbles gradually decreasing toward Stmin toward the core surface. Top right panel: ditto for 10 km. Bottom left panel: increasingly more smaller pebbles mix in between the massive inner core for 100 km core radius. Bottom right panel: for Rc = 500 km there is a more diverse mixture of pebble sizes in the core.

Fig. 6 Contributions to the total core mass for 1 and 100 km at 39 au (upper left and right respectively) and 1 and 100 km at 10 au (bottom left and right respectively) for a high resolution run of N = 105 pebble swarms. The y-axis indicates thesize range of pebbles and the x-axis the distance from the core center. All the pebbles in a certain pixel with width Δr and height Δs are counted and weighted with respect to Mt, the total core mass. The mass contribution is shown by the color bar.

4 Discussion

One of the key assumptions that goes into out model is the initial setup of the cloud, just before collapse. In reality, there will be an organic transition from a phase in which pebbles are being collected into a clump through the SI into a phase where gravity dominates and triggers the gravitational collapse that we studied in this paper. The initial setup concerns the spatial distribution of pebbles as well at the initial velocity distribution of those pebbles. The velocity dispersion of pebbles during the SI phase are unknown and not easy to obtain without currently unfeasible ultra-high resolution simulations of the SI. Processes such as gas coupling, turbulence, dust to gas feedback and random Brownian motion play an important role in determining the right dispersion profile. If turbulence is important in the cloud, small-scale clumps resulting from the SI are easily destroyed again (Johansen et al. 2011; Klahr & Schreiber 2020) where the latter authors do note that this effect vanishes likely as gas depletes in the later disk stage. The typical speed induced by turbulence is given by $v_{\text{turb}}=\sqrt{2\alpha \text{St}}{c}_{\text{s}}$ for St < 1 (Ormel & Cuzzi 2007) with α the turbulence parameter. For a standard value of α = 10⁻³ and St = 10⁻² we obtain v_turb ~ 0.5 m s⁻¹. Since the terminal velocities of the smallest pebbles are of similar order of magnitude (see Eq. (16)), we dothink turbulence may be an important effect to take into account in future follow-up studies. In general for the final collapse phase of the global SI (which we track), starting from a cloud on R_H, we believe that the terminal settling velocities are dominant due to the strong gas coupling regime that we consider. We therefore use terminal velocities to initialize the radial velocities. For the angular dispersion we use the virial dispersion profile based on the initial kinetic and potential energy of the cloud. It is not certain how realistic this assumption is, but since the gas drag rapidly diminishes the non-radial velocities, we believe that the influence on our results is small. Further investigations on velocity dispersion inside SI clumps would be important.

We use spherical symmetry with the shell approach for gravity. It is well known that small radial perturbations in the density/velocity profile of a collapsing gas cloud lead to an instability in which the density peaks at certain cloud radii (Brenner & Witelski 1998). This appears to be a consequence of forcing particles to collapse at different times for a point particle collapse, destroying the self-similarity of the free-fall solution (a simple proof is given in Appendix D). Indeed, we observe the same phenomenon in runs with narrow or single sized initial pebble size distributions. For broad size distributions we do not encounter this instability due to the much higher difference in pebble fall times. It is still unclear to us whether this instability is an artifact or an intrinsic property of the self-similar cloud collapse solution.

One physical and important effect not included in our models is that of initial cloud rotation. The presence of rotation leads to rotational support of the cloud and might prevent pebbles from collapsing to the center of mass. However, if mutual pebble collision provide an effective viscosity, the pebbles could still collapse to a single core. Results in the shearing sheet approach also show that rotation is important in answering the question if a single object forms from the collapse or a binary system (Nesvorný et al. 2010b; Robinson et al. 2020). We plan to incorporate rotation and its implications in future simulations.

The high solids-to-gas ratio in the pebble cloud can lead to entrapment of gas during collapse. Modeling compressible gas is beyond the scope of this paper. We speculate that hydrostatic effects might slow down collapse. On the other hand, the gas might rapidly escape outward before the core density becomes critical. We recognize the importance of this effect and it should be investigated further in future follow-up.

Earlier work came to the conclusion that the core of kilometer sized bodies is formed initially by mid-sized (centimeter) pebbles due to optimal energy dissipation by the gas (Wahlberg Jansson & Johansen 2017). Smaller pebbles take more time to settle, and larger pebbles do not lose enough energy through gas friction. This result is indeed expected in a virial approach in which pebbles oscillate freely through the center of mass during contraction, so that collisions dominate energy losses. The oscillation is then dampened only slowly for pebbles that hardly feel the gas, and most quickly for the pebbles with high collision rates. In our model, the virial equilibrium is not observed in any of the simulations. Instead the core forms right away from the pebbles that are most massive, before other, smaller, pebbles reach the center of mass, as the cloud becomes optically thick. The largest pebbles suffer many collisions before crossing the center of mass for an oscillation. In this way not only pebbles with the ideal damping behavior through gas friction, but also larger pebbles are incorporated into the core right away. On the other hand the smallest grains will have low terminal velocity and reach the core last, leading to a planetesimal/comet surface dominated by small particles. The gas is clearly an important aerodynamic sorting mechanism for the final core structure, particularly for a unequal initial size distribution. The effect of gas was not considered important in some earlier work due to collision timescales being much lower than stopping timescales (Nesvorný et al. 2010b; Wahlberg Jansson & Johansen 2014). However, in both these studies, the estimates of collision times was based on the assumption of randomly oriented, virial-like velocities. In our simulations, as the random velocity components get initially dampened by gas and maybe also collisions, the cloud enters into a more organized collapse where gas friction turns out to be the dominant factor. Collisions still happen, but in a realistic situation where there will be some spread in Stokes numbers, these collisions will now be dominated by large particles sweeping up smaller ones in systematic motion. Therefore, we find that gas friction remains a key ingredient in the collapse model.

5 Summary and conclusions

We summarize the core collapse model that we have developed as follows. We track the evolution of a local pebble cloud that is graviationally bound resulting from the SI. The pebbles are initially placed over a sphere with radius R_H such that the mass density of the cloud is constant everywhere. The pebbles are subject to gas drag, mutual collisions and self-gravity. The initial radial velocities of the pebbles is set to their corresponding terminal velocity. This is a reasonable assumption since any random component in the radial direction is dampened to the terminal velocity due to the efficient gas coupling of the pebbles. In the angular directions we randomize the velocities using a virial dispersion profile. Pebbles are picked from an initial MRN size distribution between St_min, St_max, in accordance with the SI.

Collisions are fragmenting or sticking, depending on the relative pebble velocities. Collisions are resolved using a Monte Carlo algorithm in which we statistically pool pebbles based on the local collision rate in the radial zone. Additionally, both collisions and advection are modeled using the representative particle approach. If the cloud’s inner one percent of mass reaches an optical depth of unity during collapse, we set the location where this happens to the accretion radius, the radius at which pebbles are considered accreted onto the optically thick core.

Self-gravity is implemented by using the shell approach: a pebble only feels the gravitational pull from pebbles closer to the center of mass as if the sum of their masses resides as a point source at the center of mass. To go forward in time, the time step for collisions and advection are compared and the time step corresponding to the best resolution is taken to proceed.

The main conclusions can be summarized as follows:

• 1.

  Fragmenting collisions in gaseous collapsing pebble clouds with St_max ~ 0.2, are rare. The critical terminal velocity v_t,* at the location where the cloud becomes optically thick lies far below the fragmentation threshold for which SI will be triggered. Growth through coagulation is negligible, except for the most massive clouds with R_c ~ 100 km and beyond.

• 2.

  Comets and planetesimals collapse toward a primordial core for which the collisions have a negligible effect on altering the initial size distribution.

• 3.

  We find that although the initial size distribution is preserved through the collapse, the order of accretion is that the aerodynamically largest pebbles form the inner core and the pebble sizes decrease toward the surface of the formed core.

• 4.

  The collapse of pebble clouds is self-similar: fall timescales are the same for the same Stokes number distribution for increasing cloud mass at given orbital distance.

• 5.

  Massive clouds are optically thick in the inner region already at t = 0, indicating that during the SI phase these regions will already be highly collisional and possibly form an inner solid core. In this phase, the size-sorting combination of gravity and drag will be less efficient, so the inner core of large objects will not be size-sorted, but will represent the initial size distribution. Only closet to the surface should we see the familiar structure of large grains inside and small grains outside.

• 6.

  Comets that have passed the star are stripped from the loosely bound small pebbles at the surface. The increasingly aerodynamically larger pebbles toward the comet core are harder to lift from the surface. We predict that the core is exposed for these comets with centimeter sized pebbles on the active surface layer. We also predict that comets that have not encountered the star still have the loosely bound millimeter-sized pebble surface layer.

Our findings support the notion that comets are primordial in nature and that comets have a systematic radial structure resulting from aerodynamic sorting of the primary building blocks. Future observations and interior analysis of a larger sample of comets are vital to test the validity of our model.

Acknowledgements

The authors thank the anonymous referee for a constructive report. We thank Marc Brouwers in particular for useful last moment comments. We thank Tom Konijn, Sjoerd van der Heijden, Michiel Min and Anders Johansen for useful discussions. R.V. acknowledges funding from the Dutch Research Council (NWO), project number ALWGO/15-01. J.D. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 714769.

Appendix A Terminal velocity collapse

Fig. A.1 Falltimes in years (y-axis) for pebble clouds composed of a fixed Stokes number (x-axis). The fall times decrease rapidly for increasing Stokes number and the numerical result agrees well with the analytical expectation given in Eq. (A.5).

Fig. A.2 Convergence test for the fiducial model for different number of representative pebbles N. Top left panel: average collisions per pebble in time. Top right panel: ratio of the kinetic to potential energy K∕|U| of the cloud in time. The results are the same for increasing N already at N = 104. Bottom left panel: CMF (settled pebbles) in time. Results are the same for all N. Bottom right panel: mass distribution function vs. the radius (Stokes number) of the pebbles of the final core.

In this appendix, we compute the collapse of a pebble cloud in the limiting case in which the radial velocities are always given by the terminalvelocity, at which gravitational force and drag force are equal in size. The terminal velocity of a pebble is then given by: $$v_{\text{t}}=-\frac{G{M}_{\text{enc}}}{r_{\text{p}}^2}{t}_{\text{s}}.$$(A.1)

For St ≪ 1, the rate of change in time of the radial position r_p of a pebble is then given by: $$\frac{\text{d}{r}_{\text{p}}}{\text{d}t}=-v_{\text{t}}.$$(A.2)

Integrating this expression $${\displaystyle {\int }_{R_0}^{r_{\text{p}}}{r}_{p^{\prime }}^2}\text{d}{r}^{\prime }=-{\displaystyle {\int }_{t=0}^t\text{d}}{t}^{\prime }G{M}_{\text{enc}}{t}_{\text{s}},$$(A.3)

with R₀ the initial release distance of the pebble and $M_{\text{enc}}~{\rho }_0{R}_0^3$ the enclosed mass at this particular release distance, gives us r_p as a function of time: $$r_{\text{p}}(t)=R_0{\left(1-4\pi {\rho }_{0}G{t}_{s}t\right)}^{\frac{1}{3}}.$$(A.4)

The time needed for a pebble to fall toward the cloud center r = 0 is: $$t_{\text{t}}=\frac{{\Omega }_0}{4\pi {\rho }_{0}G\text{St}},$$(A.5)

where we have replaced the stopping time by the Stokes number using Eq. (5). It is interesting to compare this gas-moderated collapse time with the free-fall time that would apply in the case of no friction $$t_{\text{ff}}=\sqrt{\frac{3\pi }{32G{\rho }_0}}.$$(A.6)

Fig. A.3 Same results as discussed in Fig. A.2 now for 100 km core radius. The most notable difference is a faster collision convergence due to better statistics (more collisions due to higher cloud mass).

In the presence of gas the timescale of collapse is effectively a slowed down free-fall collapse. The slow-down factor depends on the couplingstrength of pebbles to the gas (Stokes number). A comparison of the numerically determined collapse times with Eq. (A.5) is shown in Fig. A.1. The derivative of Eq. (A.5) d t_t ∕dSt ∝−1∕St² showing that the difference in fall times become increasingly larger for smaller pebble sizes. In the limit of large Stokes numbers, the gas-moderated fall time appears to become even shorter than the free fall time. However, by that time, the assumption that a particle could reach the terminal velocity is no longer valid, in effect, velocities are always limited by the free-fall velocity.

Appendix B Average collisions per pebble

To determine the importance of collisions, we estimate the total collisions a pebble undergoes before reaching the cloud center. We take an extreme case where the largest pebble i with radius s_max falls through a column of length R_H consisting of the smallest pebbles j with radius s_min. We focus on the radial collision rate fed by the terminal velocities of the pebbles. The sweep-out column (average collisions) of the pebble with s_max before settling is then given by: $${\overline{N}}_{\text{coll},i}=n_j{\sigma }_i{R}_{\text{H}},$$(B.1)

with n_j = fρ₀∕m_p,j the number density of the small pebbles representing a fraction f of the cloud mass, ${\sigma }_i=4\pi s_i^2$ the cross-section of the largest pebble. The cross-section of the smallest pebble has been omitted since s_min ≪ s_max. Writing thisexpression in more global parameters gives: $$\begin{array}{l}{\overline{N}}_{\text{coll},\text{i}}~28{\left(\frac{s_{\text{min}}}{0.01\text{cm}}\right)}^{-3}{\left(\frac{s_{\text{max}}}{0.5\text{cm}}\right)}^2{\left(\frac{r_0}{39\text{au}}\right)}^{-2}{\left(\frac{{\rho }_{•}}{1\text{g}{\text{cm}}^{\text{-3}}}\right)}^{-2/3}\\ \times {\left(\frac{M_{\star }}{M_{\odot }}\right)}^{2/3}\left(\frac{R_{\text{c}}}{1\text{km}}\right),\end{array}$$(B.2)

where we have used that the internal density is the same for the pebble and the solid core with radius R_c and the fraction of mass f the small pebbles fill in is estimated at approximately 10⁻². For the used disk radii and core radii in our model the average collisions for the largest pebble are as follows: at 39 au for 1 km, 10 km and 100 km give ${\overline{N}}_{\text{coll},i}~$28, 280, and 2800, respectively, with numerical values from the simulations at ~10, 90 and 900, respectively. For the 10 au case, 1 km, 10 km and 100 km give ${\overline{N}}_{\text{coll},i}~$53, 530, and 5300, respectively, with numerical values at ~20, 100 and 1100. The values agree reasonably well. We note that the optical depth and consequently the accretion radius lower our collision numbers.

Appendix C Convergence tests and core formation

To verify that results are independent of resolution, we run our model additionally for N = 10³ and N = 10⁵ to identify the trend of increasing the resolution. In Fig. A.2, we show that N = 10⁴ provides sufficient resolution to ensure convergence of the results. For the fiducial model, average collisions per pebble $\overline{N}$ (Fig. A.2, top left panel) are quite low, which lead to statistically noticable fluctuations. Nevertheless, the average amount of collisions per pebble is reasonably the same for different N. For the energy evolution, the graininess increases in the energy curves for lower N since a settled pebble removes a higher contribution from the initial energy of the cloud, both in potential and kinetic energy. The final mass distribution for the core is very similar for all different N (Fig. A.2, bottom right panel) indicating that collisional dynamics remains unchanged for sufficient resolution.

Additionally, we show that for the R_c = 100 km case, the convergence in resolution is also ensured in Fig. A.3. Specifically, since there are more collisions per pebble, $\overline{N}$ converges sooner in resolution due to better statistics. We focus next on the core formation procedure in our cloud collapse model.

Since no initial core is present, we quantify core formation by considering the optical depth: $${\tau }_{\text{in}}={\sum }_{i=1}^{0.01M_t}\frac{{\sum }_{j=1}^{N_z}{N}_{j,p}{\sigma }_j}{4\pi r_i},$$(C.1)

where N_j,p is the total amount of physical pebbles one representative particle with mass m_p represents, N_z the amount of representative pebbles in a zone and σ_j the cross-section of the physical pebble. 0.01M_t represents the total amount of pebbles from the inside out containing the 0.01M_t cloud mass and r_i the radial distance of the center of a zone from the cloud center r = 0. If τ_in > 1 a pebble further out cannot penetrate the corresponding layer anymore. The initial core consists of roughly one percent of the cloud mass M_t in the largest pebbles from our size range and reaches τ_in ~ 1 well before the collapse is finished. This makes sense since the largest pebbles have the shortest fall time. In the occurrence of this core formation criteria, the radius for which pebbles are considered accreted r_acc is set to the location where τ_in = 1. Simultaneously we settle the corresponding 0.01M_t and repeat the process for the next active inner region.

Appendix D Radial shell instabilities

In Newton’s shell theorem we find that pebbles in a spherically symmetric system only feel the gravitational attraction of the mass situated below them with respect to the center of the cloud r = 0, as if this mass were in the center. As a first step, we investigate if pebbles are stable to perturbations in the first place. To find out if an outer pebble situated at initial distance R_i+1 can catch up with an inner pebble with distance R_i < R_i+1, we equate their distance evolution in time given in Eq. (A.4). We assume that the system below the pebble at R_i collapses according to Eq. (A.5). We know then that the pebble at R_i also obeys this collapse time and has a distance evolution of: $$r_i(t)=R_i{\left(1-4\pi {\rho }_{i}G{t}_{s}t\right)}^{\frac{1}{3}}.$$(D.1)

The pebble at R_i+1 evolves according to: $$r_{i+1}(t)=R_{i+1}{\left(1-4\pi {\rho }_{i+1}G{t}_{s}t\right)}^{\frac{1}{3}},$$(D.2)

assuming an equal stopping time. Equating these expressions gives the time after which the two pebbles meet: $$t_{\text{t},\text{i},\text{i+1}}=\frac{R_{i+1}^3-R_i^3}{4\pi G{t}_s(R_{i+1}^3{\rho }_{i+1}-R_i^3{\rho }_i)}.$$(D.3)

Using ${\rho }_{i+1}=M_{i+1}/(4\pi R_{i+1}^3/3)$ and ${\rho }_i=M_i/(4\pi R_i^3/3)$, it becomes $$t_{\text{t},\text{i},\text{i+1}}=\frac{R_{i+1}^3-R_i^3}{3G{t}_s(M_{i+1}-M_i)}.$$(D.4)

Since we are using the enclosed mass approach, we can write M_i+1 = M_i + m_s where m_s = M_t∕N is the mass of one pebble swarm in our numerical simulation. This leads to $$t_{\text{t},\text{i},\text{i+1}}=\frac{N(R_{i+1}^3-R_i^3)}{3G{t}_s{M}_t}.$$(D.5)

Now we use the spherically uniform radial profile we use for pebbles in the cloud. Pebbles are initiated on positions $r_i={(i/N)}^{1/3}{R}_{\text{H}}$ with i the pebble number going from 1 to N. This means that a pebble put at distance $R_i=R_{\text{H}}{(i/N)}^{1/3}$ is followed by a pebble at position $R_{i+1}=R_{\text{H}}{((i+1)/N)}^{1/3}$. Filling this in recovers our assumption that pebble i + 1 should fall with the same timescale as the pebble at R_i. However, if we place the pebble i + 1 closer toward the pebble at i at $R_{i+1}=R_{\text{H}}{(i+f)}^{1/3}$, where f is a real number between zero and one and a smaller f puts the ou:ter pebble closer to R_i. Plugging this in for R_i+1 in Eq. (D.5) gives: $$t_{\text{t},\text{i},\text{i+1}}=\frac{f{\Omega }_0}{4\pi G\text{St}{\rho }_0}.$$(D.6)

This result can be interpreted as follows: the collapse timescale of a pebble with distance R_i+1 on a pebble with distance R_i < R_i+1 is always shorter than the collapse timescale of the rest of the cloud that is situated below pebble R_i. That is, if the distance separation is closer than spherically uniform between the above mentioned pebbles e.a. 0 < f < 1 (if f = 0 the collapse time is zero since they are on top of each other). In the limit of a very high f, we see from Eq. (D.6) that the collapse time becomes too high. This makes sense since they will never meet. The expression itself appears to be independent of the resolution N. This shows that the radial instabilities occur rapidly if perturbations are triggered, explaining the formation of density peaks during the collapse evolution. The pebble distribution is therefore unstable and minor changes in the initial distribution will lead to regions of high local densities. With a wide size distribution, this does not cause problems.

References

All Tables

All Figures

Fig. 1 2D sketch of the spherical cloud collapse model with the center of mass (CM) in theorigin. Pebble swarms (blue circles) are initiated over a Hill sphere with spherical velocities and positions. During evolution pebble swarms feel both gravity fg and gas drag fd. The dashed shells indicate the zones for collision evolution. As the pebbles settle, zones are being rebuilt such that the number of swarms per zone never falls below a desired value Nz.

In the text

	Fig. 2 Overview of the algorithm starting from cloud initialization at t = 0.
In the text

Fig. 3 Simulation results for three different cloud masses corresponding to core radius Rc = [1, 10, 100] km at 39 au from the central star. Top left panel: average collisions per pebble $\overline{N}$. Top right panel: ratio of kinetic to potential energy E = K/|U|, bottom left panel: CMF in time. Bottom right panel: mass weighted initial and final pebble Stokes number distribution. The initial St distribution (red curve) is the same for all three cloud masses.

In the text

	Fig. 4 Same results as shown in Fig. 3 but now for an orbital distance of 10 au for N = 104 in all three cases. The most notable changes are a significantly shorter collapse time, increasing number of collisions due to smaller Hill radii, and a gradually increasing CMF. Growth can be significant as a result of sticking collisions for 100 km core radius.
In the text

Fig. 5 2D slices in the XY plane of the final core formed from the collapse on 39 au normalized in units of Rc. The z-dimension has been flattened to one percent of the core radius. Pebbles are indicated with the circles and scale from smallest circle to largest circlewith Stmin, Stmax resp. Top left panel: core structure for Rc = 1 km. The inner core is composed of primarily Stmax pebbles gradually decreasing toward Stmin toward the core surface. Top right panel: ditto for 10 km. Bottom left panel: increasingly more smaller pebbles mix in between the massive inner core for 100 km core radius. Bottom right panel: for Rc = 500 km there is a more diverse mixture of pebble sizes in the core.

In the text

Fig. 6 Contributions to the total core mass for 1 and 100 km at 39 au (upper left and right respectively) and 1 and 100 km at 10 au (bottom left and right respectively) for a high resolution run of N = 105 pebble swarms. The y-axis indicates thesize range of pebbles and the x-axis the distance from the core center. All the pebbles in a certain pixel with width Δr and height Δs are counted and weighted with respect to Mt, the total core mass. The mass contribution is shown by the color bar.

In the text

	Fig. A.1 Falltimes in years (y-axis) for pebble clouds composed of a fixed Stokes number (x-axis). The fall times decrease rapidly for increasing Stokes number and the numerical result agrees well with the analytical expectation given in Eq. (A.5).
In the text

Fig. A.2 Convergence test for the fiducial model for different number of representative pebbles N. Top left panel: average collisions per pebble in time. Top right panel: ratio of the kinetic to potential energy K∕|U| of the cloud in time. The results are the same for increasing N already at N = 104. Bottom left panel: CMF (settled pebbles) in time. Results are the same for all N. Bottom right panel: mass distribution function vs. the radius (Stokes number) of the pebbles of the final core.

In the text

	Fig. A.3 Same results as discussed in Fig. A.2 now for 100 km core radius. The most notable difference is a faster collision convergence due to better statistics (more collisions due to higher cloud mass).
In the text