High-resolution Study of Planetesimal Formation by Gravitational Collapse of Pebble Clouds

We performed numerical simulations of the gravitational collapse of pebble clouds using the multiphysics code GIZMO (Hopkins & Raives 2015). The predominant feature of the GIZMO code is a novel Lagrangian hydrodynamics method, meshless finite mass (MFM). MFM works by initially partitioning the domain into fluid elements using a kernel function, and then continuously morphs the faces between fluid elements over time such that the mass of each fluid element remains constant. Because the fluid elements have zero mass flux, this allows them to be tracked over time, thus making MFM a Lagrangian method. The Riemann problem is then solved across the element faces with a Godunov-type scheme (HLLC method).

For the problem of gravitational collapse, Lagrangian codes are superior to Eulerian mesh codes because they conserve angular momentum, which is crucial for accurately modeling this physical process. The Lagrangian predecessor of MFM is the smoothed-particle hydrodynamics (SPH) method. The advantage of MFM over SPH is that MFM does not need an artificial viscosity term for numerical stability. This makes MFM a less diffusive method than SPH and better at preserving shock-fronts (Deng et al. 2019). Preserving the sharpness of density gradients is another important requisite of accurately modeling gravitational collapse. In our model, we definitively need to avoid viscous effects that lead to unwanted redistribution of angular momentum. The specific differences in outcome between the MFM and SPH methods for our setup are briefly discussed in Appendix C.

The initial pebble cloud is a uniform density sphere inside a Keplerian shearing box, which approximates the tidal forces in the disk on the cloud, thus in rigid rotation at the Keplerian frequency. Though this initial condition is more idealistic than using pebble clouds formed by hydro simulations of disk instabilities, as done in Nesvorný et al. (2021), it allows us to isolate the impact of the birthplace of planetesimals on their properties. The rigid rotating pebble cloud in bound rotation around the star is not only the exact condition for the concept of hill density to trigger collapse (Klahr & Schreiber 2020), it also resembles the average pebble cloud, and one can later study the effect of extra spin inside the cloud on the outcome of the collapse.

All simulations discussed in this paper are done with N = 10⁵ particles, which is the maximum feasible resolution. The resolution convergence test is discussed in Appendix D. With our model, this resolution is able to resolve tighter binaries than any previous simulations of planetesimal formation by gravitational collapse. We discuss this point further in Section 3.4. Note that our simulations are single-species, meaning we only model the pebbles and ignore the drag pressure on the pebbles from the gas in the disk during collapse. This follows Nesvorný et al. (2010, 2021), in which they found that, at high particle densities, drag forces become insignificant compared to collisions. Furthermore, as argued in Klahr & Schreiber (2020, 2021), the collapse only starts when the turbulent diffusion time of the pebble cloud is longer than its freefall; thus, during the collapse, further feedback from the gas can be neglected. Therefore the mass of a gravitationally unstable pebble cloud depends on the turbulence of the gas, which in the case of SI corresponds to masses equivalent to 10–100 km-sized solid objects.

Independent of mass, for a cloud in a circumstellar disk to collapse under its own self-gravity, the cloud's density must reach the threshold at which self-gravity overtakes the tidal forces in the disk. The initial density of our pebble cloud is set to this threshold. The traditionally used density threshold is the Hill density

where M is the mass of the central star, and a is the semimajor axis of the cloud in the disk. The Hill density is derived from the Hill sphere, i.e., the sphere connecting the two Lagrangian points at either end of the central object's Roche lobe. A Roche lobe is the teardrop-shaped region around an object in a binary system formed by the system's line of equipotential, meaning anything in that region is gravitationally bound to that object.

Note that by construction, an object's Hill sphere always extends beyond its Roche lobe. Hence, a particle can be inside an object's Hill sphere yet outside its Roche lobe, resulting in that particle not actually being gravitationally bound to the object. Because of this, the Hill density is in fact an under-approximation of the critical density for gravitational collapse against tidal shearing. We use a more suitable estimation of the critical density using the sphere of equivalent volume to the object's Roche lobe, the radius of which is derived in Paczyński (1971). This "modified" Hill density is given by

The modified Hill density is a factor of (3/2)³ greater than the traditional Hill density. Using this modified Hill density ensures that all particles in the cloud are gravitationally bound and will participate in the collapse.

It is not possible to treat all pebbles involved in the collapse as individual particles. Therefore, one treats them rather as superparticles representing a cloudlet of pebbles. But one cannot replace such a cloudlet of the combined mass of its pebbles into a solid sphere of the same material density, as one would then chiefly underestimate the collision rates of the pebble cloudlets (Nesvorný et al. 2020). And without collisions, the collapse would lead to a virialized cloud rather than a solid object (Wahlberg Jansson & Johansen 2014).

Previous N-body simulations of collapsing clouds consisted of collisional superparticles that represent subclouds of pebbles (Nesvorný et al. 2010; Nesvorny et al. 2019; Nesvorný et al. 2021). During collisions, superparticles either exchange momentum based on a collisional timescale estimator or if particles get too close to each other, then they merge to create a new spherical superparticle with a mass based on the energy of the collision and previous masses. Instead of treating collisions directly, we opted for modeling our pebble subclouds as adiabatic gas particles. The pressure of an adiabatic gas is related to its density by the polytropic equation of state:

where K is a constant related to entropy, and n is the polytropic index.

A polytropic index of n = 0.5, for instance, nicely describes the internal structure of the Earth (O'Keefe 1966). We use a slightly softer polytropic index of n = 1, which still resembles the radial density profile of Earth (see Figure 20 in Appendix A), but reflects the compressibility of a granular object and leads to a lower speed of sound, which impacts the computational time step. Our code solves gravitational equations simultaneously with the hydro equation of Equation (3), allowing pressure-supported clumps to form from the initial sparse cloud when the gas reaches hydrostatic equilibrium (HSE). The constant K can be set such that the gas stops collapsing when it reaches a density of ρ_f = 1.0 g cm⁻³, resulting in clumps with densities matching that of solid planetesimals.

We normalized our calculations by setting the gravitational constant G = 1, the central star mass M_* = 1, and the orbital radius of the cloud a = 1. This fixes the modified Hill density and therefore the initial density of our clouds, and it gives a ratio of the cloud's orbital period to freefall time of T_k /t_ff = 18. By also fixing the initial cloud radius, these simulation units allow us to use the same initial conditions for clouds with any central star mass and orbital radius. To differentiate clouds at varying orbital radii around a star of a given mass, we define the ratio of the initial cloud radius R₀ to the collapsed equilibrium radius R_f , defined for the case that all of the pebbles would end up in a single body, as the collapse ratio η_c

In physical units, the final collapsed density of clumps is set to ρ_f = 1.0 g cm⁻³, and the initial density of the cloud is given by ${\rho }_{{H}^{* }}$, which depends only on the central star mass M_* and the orbital radius of the cloud a. This means that from cloud to clump, the density will increase by a factor of ${\left(\tfrac{{\rho }_{{H}^{* }}}{{\rho }_{f}}\right)}^{-1}$ and the radius will decrease by a factor of ${\left(\tfrac{{\rho }_{{H}^{* }}}{{\rho }_{f}}\right)}^{(1/3)}$ (assuming that the density is constant over the volume). Therefore, where the cloud collapse is occurring in the parameter space of orbital radius and central star mass can be entirely described by the collapse ratio η_c .

The HSE density of clumps in our calculations is set by K, so to use η_c as a variable parameter, a relation must be found between η_c and K. Using the Lane-Emden equation to solve for the HSE value of K for a clump with ρ_f = 1.0 g cm⁻³ gives the following equation for K:

where G is the gravitational constant, R₀ is the initial cloud radius, ρ₀ is the initial cloud density, ξ₁(n) is a dimensionless radius variable, and $\tfrac{{\rho }_{c}}{\bar{\rho }}(n)$ is the ratio of the central to mean density of the final clump, for which we choose ρ_c = ρ_f . The terms ξ₁(n) and $\tfrac{{\rho }_{c}}{\bar{\rho }}(n)$ come from the Lane-Emden equation solution. The full derivation for Equation (5) is given in Appendix A.

We assume that all of our initial pebble clouds are in bound rotation around the central star, which follows from the ansatz of Hill density, and thus they will all have the same centrifugal radius r_c to Hill sphere radius ratio. The centrifugal radius for our pebble cloud is given as

as we derive in Appendix E. This means that if η_c ≥ 0.13, then a single object can form from the collapse, which may be the case inside 0.02 au. For large distances, r_c will become much larger than R_f , and thus binaries will be born with larger separation, and multiple fragmentation becomes more likely. This is the scientific question of this paper: the properties of binary and higher-order planetesimal systems as a function of their distance to the Sun.

Ideally, for the correct pressure–density relation of solid planetesimals, the Tillotson equation of state should be used in combination with solving the equations for elasticity and plasticity of the body (Schäfer et al. 2016). However, this is beyond the scope of the current paper, in which we want to study the phase transition from collision dominated pebble cloud regime to a quasi-incompressible body.

The interaction between two mass elements in the GIZMO code is described via the acceleration according to the local pressure gradient

with pressure P = K ρ^γ , where $\gamma =1+\tfrac{1}{n}$ and sound speed ${c}^{2}=\tfrac{\partial P}{\partial \rho }$ . For ρ = ρ_solid, i.e., when the body should be in HSE, one can show that c is on the order of the escape velocity of the formed planetesimal. In the dilute phase, when ρ < ρ_solid, c² = γ K ρ^γ−1 is much lower than the typical relative velocities v of the pebble clumps, as long as γ > 1. Thus, shocks will occur, which mimic collisions in our model.

Deviating from the continuity approach, we can consider Equation (7) as describing a repulsive force between two particles of mass m and volume proportional to the cube of separation r for a "subsonic" collision c² > v². We multiply both sides by ρ and integrate over the volume V ∝ r³ and get a formal equation for the acceleration of a superparticle with mass m from the interaction with a second superparticle of the same mass,

More rigorously, we could use the SPH approach of treating the pressure gradient in GIZMO in which the acceleration is proportional to pressure times the gradient of the kernel function W (Springel 2005), which also leads to Equation (8). In a two-body collision, such a repulsive potential of shape ${c}^{2}\mathrm{ln}r$ would perfectly model the elastic momentum transfer from one particle to the other and well approximate the collision of two hard spheres for sufficiently long time scales. Thus, a subsonic collision between two of our superparticles is already sufficiently modeled. With our polytropic equation of state, which conserves entropy, we then find that subsonic collisions are perfectly elastic.

Nesvorny et al. (2019) emphasized the importance of inelastic bounces among superparticles representing clouds of pebbles to efficiently cool their rms velocities, which are powered by virialization during collapse (Wahlberg Jansson & Johansen 2014). At first glance, a polytropic equation of state implies perfect reversibility of any compression. Thus, there should be no damping of kinetic energy. Fortunately, the GIZMO code solves the Riemann problem for interacting particles, which only conserves entropy if the involved velocity difference is subsonic. For the densities at the start of collapse, the sound speed c can be compared to a typical collision velocity, for which we choose the freefall or escape velocity of a body with mass m and density ρ. We calibrate the pressure, or more precisely K, for the formed planetesimal at size R_f and solid density ρ_f (see Equation 5):

Thus, we can express the entropy via the escape velocity

and directly express the speed of sound in terms of the escape velocity

Here we see that for the body at solid density ρ_f , the speed of sound is of order unity, identical to the escape velocity. For lower densities, i.e., at the start of the collapse, the speed of sound is given as

which can be compared to the freefall velocity during collapse:

implying that the Mach number of a typical collision scales as:

So if our initial density is ρ = 10⁻⁹, typical collisions are highly supersonic (M = 1000), and thus dissipate enough energy to form planetesimals, all while momentum, and more importantly angular momentum, is perfectly conserved. A collision with an impact parameter larger than zero is then an off-center flyby of dilute superparticles and happens typically at supersonic velocities c² < v². Here the Riemann solver automatically replaces the gas pressure by the ram pressure:

From that, we can calculate the loss of kinetic energy in a collision as

in which ${t}_{\mathrm{coll}}=\tfrac{r}{v}$ is a representative collisional time, giving a physical meaning to the decay of kinetic energy.

In classical collision schemes (Nesvorný et al. 2020) for superparticles (SPs), one starts with treating a system of real particles (RPs). Using a function of the RPs' size distribution dN and the cross section of the individual pebbles $\sigma ^{\prime}$, one can then calibrate a new cross section σ for the SPs, which each contain n RPs that receive the same loss of kinetic energy per time. Nesvorný et al. (2020) found that $\sigma =\sigma ^{\prime} n$ is a good approximation for the case that all SPs have the same mass:

as is the case in our simulations. Thus, we can write our prescription as:

to be identical to the classical collision scheme as long as $\tfrac{\sigma ^{\prime} n}{{r}^{2}}$ is of order unity. We mimic the SP collision scheme without defining the actual number n of RPs in our SPs, which would be given implicitly by their radius a₀ and the current cross section for the SPs $\sigma \propto {a}_{0}^{2}$ as $n=\tfrac{{r}^{2}}{{a}_{0}^{2}}$. So, as soon as "solid" planetesimals of 100 km in diameter form and r ≈ a₀, then we get n = 1.

For a proper treatment of the pebble cloud collapse, the ratio between collision time and freefall time must be correctly captured. The collision time for the initial pebble cloud at virial velocity is given as

which produces the same result for different pebble sizes a₀ as long as ${{na}}_{0}^{2}$ is kept constant. We can also express the collision time as a function of η_c :

and we find that for a typical η_c = 10⁻³, the entire expression $\tfrac{{a}_{0}}{{R}_{0}}{\eta }_{c}^{-2}$ equals unity for a cloud mass equivalent to a 100 km planetesimal consisting of 10 cm pebbles. The collision time is then

For our collision model using the shock treatment of the Riemann solver, we got ${{na}}_{0}^{2}={r}^{2}$, which produces the correct collision time at the onset of collapse r = R₀:

During the collapse, the effective cross section in our model shrinks continuously to the physical cross section, when the final planetesimal forms. This shrinking of n can be regarded as coagulation during the collapse, which otherwise is subject to modeling, because the cross section should not become larger than what is given by the square of SP separation.

The effective pebble size and cross sections during collapse scale as:

and for η_c = 1 we have n = 1 and a₀ = r₀ with r₀ the size of one element at solid density. That size depends on the number of particles N we simulate and the compressed size of the pebble cloud as ${r}_{0}={R}_{f}{N}^{-\tfrac{1}{3}}$. Also mass has to be conserved; thus, ${{na}}_{0}^{3}={r}_{0}^{3}$. We combine this into

for R_f = 50 km and for η_c = 10⁻², that is a₀ = 10 cm, and for η_c = 10⁻³ that is a₀ = 0.1 cm. So at 30 au at the onset of collapse, the pebbles are a little too small, but during collapse they quickly "grow." Using too few particles would then lead to an even shorter collision time, and using too many particles would underestimate collision, so one has to choose these numbers carefully. One can imagine compensating for this inside the Riemann solver, with a modification of the RAM pressure, but this shall be left to future work.

A more detailed collisional model would also have to include the size distribution of pebbles constituting the pebble cloud and even allow for pebble clouds of different mass. Also, the coefficient of restitution C cannot be meddled with in our setup. As we solve no explicit energy equation, we effectively have very inelastic collisions C ≪ 1. To summarize, our collision model neither defines a coefficient of restitution C nor the typical pebble size a₀. Yet by construction it reproduces the correct collision timescale for the regime of interest, i.e., the collision time to be shorter than the freefall time. If one would wish to simulate other cases, like much longer (or shorter) collision times, one would have to modify the Riemann solver, such that the RAM pressure would be artificially lowered (or increased). Regardless, this is beyond the scope of the present paper.

We performed tests of collapse simulations for various resolutions and found the dissipation of kinetic energy to be virtually independent for the range of resolutions we picked (see Figure 23 in Appendix D). This is a necessary condition for a superparticle collision model (Nesvorny et al. 2019), which is met by our pressure and Riemann solver method for handling collisions. This approach intrinsically handles the collisions of pebble subclouds while also allowing the planetesimals to develop a physical size as they form, unlike the collisional superparticle method.

Despite the ad hoc ansatz of using pressure as a proxy for collisions, it has overwhelming benefits for the performance of numerical simulations. It treats collisions in a qualitatively correct way, dissipates kinetic energy lost by inelastic collisions, and it also provides a smooth transition from compressible to a quasi-incompressible planetesimal. The transition crosses up to 10 orders of magnitude in density in our simulations, and still the model provides a smooth adaptive resolution for the entire process.

To identify clumps, we use the Scikit-learn (Pedregosa et al. 2011) implementation of the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm (Ester et al. 1996). DBSCAN is preferable to other clustering algorithms such as K-Means because it does not require the number of clusters to be specified. This allows the finder to be general enough to handle analyzing collapses that produce variable amounts of clumps. DBSCAN takes two input parameters, a distance and a minimum number of points, that act as a density threshold. The specifics of how DBSCAN works and how we use it are explained in detail in Appendix B. Before passing data to DBSCAN, we take an additional step and remove particles below a certain density threshold ${\rho }_{\min }=0.1 \% {\eta }_{c}^{-3}{\rho }_{{H}^{* }}$. Particles below this density are categorically not in clumps, and filtering them out improves both the efficiency and accuracy of the DBSCAN algorithm. In short, the clump finder routine we have built defines a clump as meeting the following three criteria:

• 1.  
  All points in a clump are > = 0.1%ρ_f
• 2.  
  Clump masses are > = 0.1% M_cloud
• 3.  
  All points in a clump are within a distance R_f from at least one other clump point.

In effect, this defines the smallest clumps that we will identify to be 0.1R_f . These clumps are still represented by up to 100 individual GIZMO particles, yet they will not contribute significantly to the cumulative mass distribution (CMD) of objects, as we will show a posteriori.

Once the individual clumps are identified, we determine the binaries and higher-order systems in much the same way as Nesvorný et al. (2021). We perform a simple double loop over all of the clumps, starting with the most massive. A binary is confirmed if the two clumps have a negative binding energy and are within a Hill radius of each other. Then, an additional loop is performed over the binaries to find satellites of binary systems.

After the initial collapse phase when the particles become very dense, the evolution of the formed clumps become significantly more computationally expensive. The time step is limited by the sound crossing time in our formed planetesimals, which is independent of the heliocentric distance a of the particular system, yet the system evolves on the freefall and dynamical time scales at the given distance to the Sun, which both scale with t_dyn ∝ a^1.5. Therefore a simulation at 30 au is more than 460 times more costly than one at 0.5 au. In order to compare our results, we stopped all simulations at 1.43 freefall times. At this time, between 80 and 95% of the pebbles are typically accreted into planetesimals (See Figure 8). Ending the study here is justifiable in the sense that the scope of this study is only meant to characterize the properties of the "initial" planetesimal population that forms from the gravitational collapse of the pebble cloud. Of course one would also be interested in the longer-term dynamical evolution of the planetesimal clusters, but this has to be done by other numerical means, like with a sink-cell approach. Also, after a few freefall times, the planetesimal cluster can no longer be considered as isolated from the protoplanetary disk environment, for instance other planetesimals. Therefore, to study the long-term evolution of the planetesimals, the external disk environment must be modeled and accounted for in dedicated N-body simulations. Pushing our simulations along would only result in a further dispersal of the initial population and perhaps some mergers. The future goals of this work are to take the initial planetesimal populations we formed in these simulations and insert them back into a disk environment to more accurately study their evolution and survival rate.

We simulated cloud collapses across the solar system, from 0.5–30 au, in order to determine how the birthplace of planetesimals via the collapse ratio η_c affects their properties and multiplicity. The orbital radius proved to have a profound influence on the dynamics of the collapse and the subsequent properties of the formed planetesimals. Figure 2 provides a general look at the dynamics of a cloud collapse in a circumstellar disk. Beginning with a uniform sphere, the cloud begins to collapse into an ellipsoid shape due to the tidal forces squishing the cloud in the z-direction and radial direction. As the ellipsoid continues collapsing, first a disk and then a high-density bar forms. Qualitatively, these are the same stages of collapse that were also reported on in the work by Nesvorný et al. (2021). The tidal forces rotate the bar into a spiral, and several high-density subregions within the bar collide along with premature clumps. The spiral does one more quarter rotation, which marks the end of the dynamical collapse. Clumps that have formed are left to interact with the other clumps and accrete the residual pebbles from the collapse on secular time scales.

Zoom In Zoom Out Reset image size

Though the overall collapse proceeds similarly in both the 3 and 30 au case, there is one key difference between them, as seen in Figure 2. Throughout the entire collapse, there is a significantly higher degree of fragmentation in the 30 au cloud. This manifests as both an increase in the number of forming clumps and a decrease in their average size. Similarly, the lesser amount of fragmentation results in fewer but larger clumps in the 3 au collapse. Figure 3 shows the initial clump formation states of the entire suite of collapse simulations at different heliocentric distances, with formed planetesimals indicated by faint yellow circles. This figure demonstrates the progression of fragmentation with distance from the Sun. The rest of this section analyzes the properties of the formed planetesimals as a function of disk birthplace.

Zoom In Zoom Out Reset image size

Once planetesimals form from a cloud collapse, they mark the starting point for embryo formation. Unless some of the forming planetesimals are large enough to go straight into the pebble accretion phase (Ormel & Klahr 2010; Johansen and Lambrechts 2017; Ormel 2017; Drazkowska et al. 2022), which would be typically several 100 km, they first have to continue to grow via mutual nondestructive collisions, in an oligarchic phase. In this phase, the size distribution of planetesimals stays narrow in shape, but shifts slowly toward larger sizes. Eventually one reaches a runaway phase (Morbidelli et al. 2009; Kobayashi et al. 2016), at which individual planetesimals quickly outgrow their peers. Because of this, constraining their initial mass and size distribution is essential for studying planet formation. In Klahr & Schreiber (2020), we derived a typical cloud mass that can undergo a gravitational collapse, and in the present project, we want to learn what the overall efficiency and final outcome in mass function of that collapse can be.

The size and mass distribution of the planetesimals we formed and their dependence on orbital radius are depicted in Figures 4 and 5. The size of each clump is determined by using its mass and applying a final uniform density of ρ_f = 1 g cm⁻³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figure 4 we can see that regardless of disk birthplace, the largest ≈3 planetesimals are always greater than 0.5 D_eq, the equivalent diameter of the collapsed pebble cloud. Equivalently, from Figure 6 we see that these same handful of large planetesimals contains a majority of the initial cloud mass.

Zoom In Zoom Out Reset image size

For many years, observations have shown that almost every minor solar system object family, from the asteroid belt to Kuiper Belt, has the same characteristic size of ≈100 km. Because these objects are planetesimal remnants, this implies that planetesimals have a preferred formation size regardless of their distance from the Sun. Our results in combination with Klahr & Schreiber (2020) reproduce this phenomena, strongly suggesting that gravitational collapse is the predominant route to planetesimal formation.

This will then give us the size of the original pebble clouds that form the planetesimal population. Morbidelli et al. (2009) found that most planetesimals are born big, with sizes in the range 100–1000 km. This size range was further confirmed later by Delbo' et al. (2017), and using those results we determined our equivalent diameter to be D_eq = 200 km.

The size and mass of the planetesimals relative to the entire population are shown more clearly in Figure 5. Bottom to top represents the growing distance from the Sun. In this direction, generally there is a small decrease in the mass and size of the largest planetesimal. The largest planetesimals also tend to become more equal in size as the orbital radius increases. In all of the simulations, there seems to be a gap between the large and small planetesimals occurring in the range 10⁻²–10⁻¹ M_tot. There is also a steady increase in the number of small planetesimals (<0.5 D_eq) as orbital radius grows, with a significant difference for <5 au than for ≥5 au.

This distinction for the size distribution for the inner and outer solar system can be seen more clearly in Figure 7. The objects formed in the five collapse simulations at 0.5–4 au and the five at 5–30 au, respectively, were collected into two groups and their size distribution was plotted. In both regions, the largest planetesimals are born big. In the inner disk, however, they are born slightly "bigger." In the inner region, the largest planetesimals are 95% of D_eq, whereas in the outer part they are more like 75% of D_eq. Both regions have the same distribution shape for large objects, but the outer region has a more pronounced tail in smaller planetesimals.

Zoom In Zoom Out Reset image size

Delbo' et al. (2017) identified a primordial family of asteroid fragments that is also very large and old. By eliminating this and other families of fragments in the same semimajor axis region, they could identify a number of asteroids that they claimed cannot be included in any collisional family that had formed in locally in the main belt. Their logic consequence is that these asteroids that cannot be family members are some of the original planetesimals. From these asteroids, the authors produced a size distribution of the original planetesimal population. In the size distribution of KBOs and asteroids, there is a turning point at D ≈ 100 km, and in our simulations, we see the same turning point at D = 0.5 D_eq. This implies that D_eq ≈ 200 km.

Figure 6 shows the cumulative mass of the planetesimals formed at different orbital radii in units of the initial cloud mass. There is a small but definite spread in the mass of the largest planetesimal, with the inner to outer disk corresponding to the largest and smallest. For all but the 0.5 au case, the masses of the next largest planetesimal all follow roughly the same slope. Then, for distances below 5 au, the distribution flattens to lower masses, with a fairly large jump in the mass of the next clump. Above 5 au, all of the populations show the same mass distribution of clumps below 10⁻¹ M_tot. Between 5–10 au, there seems to be a distinct regime change in the initial mass function (IMF) of planetesimals. This makes sense considering there is a significant difference in the characteristics of objects in the solar system that are above and below 5 au.

Note that in our analysis, we include all formed planetesimals, bound and unbound. This means that if the simulations were continued, planetesimal mergers/accretion could change the overall distributions shown. However, to correctly constrain the effect of future mergers on the distribution evolution, the disk environment must be considered.

The clump properties analyzed by the previous plots in this section are compiled in Figure 8 where the following clump properties are plotted against orbital radius: number, total mass formation efficiency, and the maximum, minimum, and mean of the masses and sizes. This plot provides an overview of the trends we have identified for planetesimal formation properties depending on disk birthplace. One additional insight we can make from this plot is that for all orbital radii studied, the amount of mass from the initial cloud that ends up in clumps is >75%. This confirms that gravitational collapse is a very efficient process in converting pebbles to planetesimals.

We analyzed a suite of 12 simulations that differ only in the random seed used to generate the initial positions of the particles in a uniform spherical distribution. This was done to measure the degree of stochasticity in a collapse and to see how much it affected the formation properties of planetesimals. The 12 simulations discussed in this section are at an orbital radius of 1 au. This distance was picked for the statistical study because it is the computationally cheapest collapse that still produces multiple clump systems.

The analysis covered in this section matches exactly the one in Section 3.1. This is done to serve as a benchmark for the consistency of the properties we saw in the previous set of plots. We can use the following results to determine which characteristics of the size and mass distribution of planetesimals formed across a disk are fundamental versus consequential of random perturbations in initial conditions.

Zoom In Zoom Out Reset image size

Figure 9 gives a qualitative overview of the degree to which the random seed affects the final morphology of the collapse. Each of the snapshots are taken of the orbital (x–y) plane at 1.1 t_ff when the clumps in each simulation first form. At first glance, though the orbital dynamics of the planetesimals vary (i.e., spatial distribution of the clumps), the number formed seems to be consistent at 3 ± 1.

Zoom In Zoom Out Reset image size

Figures 10 and 11 show alternative views of the clump mass and size distribution in each of the 12 seed simulations. Most of the collapses produced three to four planetesimals (67%), with a few forming five to six (25%), and one forming 12 (8%). All simulations produced two to three large planetesimals with sizes ⪆0.5 D_eq. Though the number of clumps formed seems to be fairly stochastic, the mass of the largest clump in each of the 12 simulations is roughly consistent with an average value of 60% M_tot and a standard deviation of 9.6% M_tot. The amount the largest clump mass deviates seems to depend on how comparable in mass the second largest planetesimal is. Because the mass of the largest clump is consistent for collapses with different seeds, the inverse relation we found between largest clump mass formation efficiency and orbital radius is valid even considering the stochastic nature of gravitational collapse.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 shows the combined size distribution of all of the clumps formed at 1 au for the 12 different seeds, as well as the two size regimes and power law for the original planetesimal population of the asteroid belt derived in Delbo' et al. (2017). The size regime above 0.5 D_eq again confirms the same steepness as the Delbo' et al. (2017) data, confirming this characteristic for collapse across the disk. The smaller planetesimals at 1 au have a flatter distribution than the original planetesimal population, confirming the conclusion made from Figure 7 that the slope of the small planetesimal (<0.5 D_eq) size distribution steepens with increasing orbital radius.

Zoom In Zoom Out Reset image size

Figure 13 shows the CMD of the planetesimals. In comparison to the CMD for planetesimals at different orbital radii (see Figure 6), the same two regions of growth can be seen for masses above and below 10⁻¹ M_tot. For clumps above 10⁻¹ M_tot, all seeds have the same shape with a slight spread in masses. For masses below 10⁻¹ M_tot, there are three different patterns depending on how many clumps have formed. The CMDs with a smaller amount of clumps have a fairly steep slope with the smallest planetesimal being relatively large. Intermediate-number CMDs have a flatter shape for the smallest planetesimals. And most notably, the CMD of the simulation that produces 12 clumps starts to resemble the same distinctive shape of the CMD seen for orbital distances at 10 au and above. The fact that even at 1 au, when enough small planetesimals form, the CMD matches the CMD for >10 au confirms suspicions of a universal IMF for planetesimals.

Zoom In Zoom Out Reset image size

The compilation of clump properties formed at 1 au with 12 different seeds is shown in Figure 14. The clump mass formation efficiency remains high regardless of seed, with all 12 simulations converting ≥88% of the initial cloud mass into planetesimals. In addition to the mass of the largest planetesimal being roughly constant, the average planetesimal mass is also consistent despite the stochasticity of the collapse. This spread in mass is seen more clearly in the average planetesimal sizes, where the average cloud will produce planetesimals of size 102 km on average with a standard deviation of 52 km. This is with D_eq = 200 km. The average size of planetesimals decreases as the total number of them increases. This is consistent with the previous conclusion that a higher degree of fragmentation leading to more individual clumps forming will result in more mass being distributed to smaller planetesimals. This is the same trend we saw with the properties of clumps at different orbital radii, where the farther out a collapse occurs, more planetesimals form that are on average smaller in size.

Zoom In Zoom Out Reset image size

To summarize, the stochastic nature of gravitational collapse does result in some variations in the mass distribution of planetesimals. However, the global properties of the planetesimal populations are ubiquitous despite small changes in the initial conditions of the cloud due to the use of different random seeds. This means that the overall trends and properties of planetesimal populations we identified in Section 3.1 are valid as inherent characteristics of collapses at the given orbital radius rather than random consequences of the initial setup.

In the previous section, we did a qualitative comparison of our simulation data with the observational data in Delbo' et al. (2017) by defining a set of various power laws to characterize the distribution. This type of analysis is done not only on observational data, but also in the analysis of the mass distribution of self-gravitating pebble clumps generated in simulations of planetesimal formation (Simon et al. 2017; Liu et al. 2020). In this section, we try a different approach to describe the size distribution of asteroids. We start from the idea that planetesimals form with a preferred mass, containing a certain fraction of the initial pebble cloud and that there is a Gaussian variation around this value. For this purpose, we first define the cumulative asteroid mass per mass bin:

In the case that the mass distribution of asteroids is a Gaussian distribution of asteroid masses around median mass M₀ with a width of σ_M :

then the corresponding cumulative function is the error function. Once the cumulative mass function for the asteroid data by Delbo et al. (2019) is calculated (see Figure 15), one can identify the most likely mass, M₀, where the cumulative function reaches 50% of the total mass in objects. In this fit, we removed the two largest objects in this set (Vesta with a diameter of 521 km and Flora with 220 km ⁶ ), which we suppose to be outliers of the distribution. This exclusion is mainly motivated by the fit of the Gaussian distribution, and thus is justified a posteriori. In the case of the formation of Flora and Vesta, we argue that they were not primordial planetesimals, formed directly by pebble collapse, but that they formed by collisional evolution from smaller (≈100 km) precursors. In Voelkel et al. (2021a), we show how planetesimals grow from initial 100 km sizes via mutual collisions, before at some size, pebble accretion (Ormel & Klahr 2010) takes over. Depending on where Vesta was formed and how big pebbles actually are, its size may still be smaller than expected for efficient pebble accretion. Furthermore, if Vesta would have reached a size large enough for pebble accretion while the pebble flux in the disks still lasted, it would have quickly grown to planetary embryo sizes. Furthermore, in the case of Flora, it is also possible it formed from the collapse of a higher-mass pebble cloud.

Zoom In Zoom Out Reset image size

Thus by excluding Vesta and Flora from the data set, we find the peak in the mass distribution at M₀ = 125 km. The width σ_M can likewise be found where the cumulative distribution reaches a value of $84 \% \approx 0.5({erf}1/(\sqrt{2}{\sigma }_{M})+1)$. Using the data from Delbo et al. (2019), we find a width of σ = 60% of the mean mass (see Figure 15). Overplotting the error function of the Delbo et al. (2019) data, we see a nice fit. From this, we see that in fact the primordial asteroids fall on a narrow mass distribution around 125 km-sized objects. When translating the error function in mass to the classical number per diameter plot (see 15) we also find the smaller objects to be well described in the wings of the Gaussian.

Additionally, we included our 1 au simulation data (blue plus signs) and find they also fit the mass Gaussian nicely. This allows us to derive an initial pebble cloud mass corresponding to a 152 km asteroid. Thus, the entire mass distribution of primordial asteroids can be explained by the stochasticity of the collapse from a single pebble cloud mass. Any additional masses of unstable pebble clouds as found in large-scale simulations of the SI (Liu et al. 2020) would widen and alter the mass distribution of planetesimals, albeit to what degree would have to be constrained elsewhere.

So the 152 km equivalent size of the pebble cloud can explain all primordial asteroids, except for Flora and Vesta. If they would also have been the product of a single pebble cloud collapse, the mass of this cloud would have too strong of outliers, and it would be hard to explain why this massive cloud would not also have produced smaller companion asteroids showing up in the distribution.

By interpreting our numerical data from 12 collapse simulations forming 56 planetesimals, we can interpolate that the 26 asteroids in the Delbo et al. (2019) data are the outcome of roughly six pebble cloud collapses all of similar size. In other words, each pebble cloud in the asteroid belt may have spawned on average four asteroids of the primordial population.

So far we have characterized the properties of individual planetesimal clusters stemming from a single pebble cloud. In this section, we transition into discussing the properties of the binary and satellite planetesimal population. All of the binary analysis is compiled into Figure 16, in which the following properties are considered: size ratio, binary semimajor axis, angular momentum, eccentricity, and inclination. These plots were chosen specifically to compare our binary properties to those formed by the A12 simulations in Nesvorný et al. (2021). Also shown by the plus signs in Figure 16 are the observed properties of known binary asteroids and KBOs, cataloged by Johnston (2019). The color of all of the markers represents the heliocentric semimajor axis of the system in astronomical unit.

Zoom In Zoom Out Reset image size

Of the binaries analyzed in Figure 16, the opaque circles are the tight binaries that formed with periods short enough to have their survival substantiated within the time frame of our analysis. Due to the high resolution of our study, the time span of our simulations is too short to be certain of the evolutionary survival of the wide binaries and satellites that have formed, shown by transparent squares and circles, respectively. These wide binaries and satellites that formed in collapses beyond 4 au are distinguished in the plot because their surviving properties cannot be constrained within the time frame of our simulations. This is because the long-term collisional survival of binaries in the outer solar system decreases with smaller object sizes and larger initial separations (Nesvorný et al. 2021). Therefore, the properties shown in Figure 16 reflect tight binary properties and only the initial wide binary properties.

The size ratio of binaries we produced favors equally size binaries with R₂/R₁ > 0.5 matching the observations of KBOs that also contain a majority of equally sized binaries. The outer disk binaries (>4 au) match the tightest observed KBO binaries, occupying the same regions in the parameter space of scaled binary semimajor axis versus size ratio. The observed outer disk KBOs with scaled semimajor axes around 50 do match the wide binaries we have formed. A few satellites even venture into the very wide region of the plot. The inner disk binaries on the other hand have roughly the same to slightly lower scaled semimajor axis and have a much higher size ratio than the observed binaries in the asteroid belt. This suggests that the survival rate of binaries in the asteroid belt is much lower than in the Kuiper Belt. It seems that regardless of location in the disk, binaries favorably form with equally size ratios, and any deviation from that comes from its evolutionary lifetime in the disk.

In the angular momentum plot, our binaries only exist in regions where actual observed KBOs are. The simulated outer disk binaries follow the same relation as the observations where angular momentum increases with size ratio, as well as multiplicity increasing with both. While the outer disk binaries agree with the outer disk KBOs, the inner disk binaries have a much higher angular momentum than their counterparts in the asteroid belt. This again suggests that outer disk binaries do not have to undergo many evolutionary changes over the course of their lifetime, while inner disk binaries are more prone to property-altering processes. It appears that inner disk binaries lose angular momentum during disk evolution more efficiently, which is plausible due to many possible mechanisms like gas drag, radiation pressure, or secular gravitational effects, all of which are more severe in the inner part of the solar nebula (Lyra et al. 2021).

On the eccentricity plot, there are two lines denoting the pericenter distances β = a_b (1 − e_b )/(R₁ + R₂) = 35 (black, dotted) and 2 (black, dashed). These represent the binary separation resolving limits of our model (β = 2) and the previous model (β = 35) used by Nesvorný et al. (2021). This means our models can detect and follow binaries that are 18 times closer than previously possible, essentially reaching the maximum resolving limit because we can detect contact binaries. The only caveat is that we need higher resolution and a better pressure prescription for solids to track the evolution of contact binaries. Nevertheless, we see a direct relation between increasing eccentricity and separation. As orbital radius increases, eccentricity increases along with binary separation. These same trends can be seen in the observations of KBOs. However, it seems that the eccentricity of close binaries decreases over the course of their lifetime in the solar system, and we do see this decrease in binary eccentricities begin to happen during the time frame of our study. The mechanism for this can be pebble accretion or interactions with satellites and/or other planetesimals passing by.

The inclination of every tight binary formed hardly deviates from the initial inclination of the collapsing cloud, with all deviating by less than 1°. Despite the relatively small range of inclinations, there is a clear trend that inclination increases with distance to the Sun, related to the fact that more equally massed objects provide more degrees of freedom to distribute angular momentum. We also see some satellites reaching inclinations as high as 6°. Nesvorný et al. (2021) studied the collapse of clouds with a distribution of spin orientations from SI simulations using an overall spread of inclinations over 360°. Still, the deviation of inclination between cloud and binary was at most in the range of Δi = 10°–20°, suggesting that either wider binaries tend to vary more in inclination, or that the initial spin of the pebble cloud has a profound impact on the forming binary, which we cannot test in the present paper.

An important aspect of planetesimal formation by gravitational collapse is the redistribution of angular momentum to many objects through multiple fragmentation. This dispersion of angular momentum is what allows the tight binaries we observe in our solar system to form. Figure 17 shows the proximity factor ζ = 2a_b (1 − e_b )/(R₁ + R₂) which indicates the tightness of the binary, against the binary's distance from the central star. The minimum proximity factor is ζ = 1, which occurs for contact binaries. The maximum proximity factor ζ₀ is dependent on the heliocentric semimajor axis and is denoted in Figure 17 by the dashed black line. This value represents the proximity factor for a binary that would form if all of the angular momentum and mass from the initial cloud went into a single equal-mass binary (see Appendix E). From Figure 17, all of the simulated binaries and almost all of the observed binaries have ζ < ζ₀. This suggests that the majority of KBO binaries originally formed alongside other binaries and planetesimals in a multiply fragmenting collapsing cloud. In both the simulated and observed binaries, the binaries tend to get wider as the maximum proximity factor ζ₀ increases. However, most binaries are well below the ζ₀ limit, and the farther out binaries are relatively farther away from the maximum proximity factor, indicating a higher degree of fragmentation in clouds in the outer disk.

Zoom In Zoom Out Reset image size

Because we fixed the initial modified Hill density of the cloud, our results can be applied to disks with different central star masses at a certain orbital radius.

The conversion from the collapse ratio η_c to a physical location in the disk can be seen in Figure 18. Note that there are multiple modified Hill densities for a given (M_*, a), meaning simulations with the same collapse ratio are applicable at multiple combinations of central star mass and orbital radius.

Zoom In Zoom Out Reset image size

This means that to apply our model to exoplanetesimal formation in extrasolar disks, the only parameter to consider is the mass of the central star. To do this, we can use Figure 19 to convert our results for a 1 M_⊙ star to another star mass. For instance, consider our 30 au cloud collapse around the Sun. This corresponds to an η_c value of roughly 7 × 10⁻⁴. This means that a cloud collapse at 30 au in our solar system corresponds to a collapse at ≈15 au in a 0.1 M_⊙ disk and ≈60 au for a 10 M_⊙ disk.

Zoom In Zoom Out Reset image size