EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 683 (March 2024) A&A, 683 (2024) A38 Full HTML
Highlight
Open Access
Issue
A&A
Volume 683, March 2024
Article Number A38
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202347742
Published online 04 March 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Planetesimals are believed to form by the gentle gravitational collapse of aerodynamically concentrated clumps of pebbles (pebble clouds). Pebbles are millimetre- to decimetre-sized aggregates that have grown by collisions of micrometre-sized dust and ices grains in the gaseous environment of protoplanetary discs. Laboratory and numerical studies of collisional dust growth have shown that bouncing, erosion, fragmentation, and radial drift limit the maximum sizes to which pebbles can grow to the millimetre- to centimetre-size range (Blum & Wurm 2008; Güttler et al. 2010; Zsom et al. 2010; Krijt et al. 2015; Blum 2018; Lorek et al. 2018). The streaming instability operates well for these pebble sizes, a process that has been demonstrated numerically to be efficient in concentrating large amounts of pebbles into filament-like structures that subsequently fragment under their own gravity to form pebble-pile planetesimals (Youdin & Goodman 2005; Johansen et al. 2007, 2014; Wahlberg Jansson & Johansen 2014).

The minor bodies in the Solar System, such as asteroids, comets, and Kuiper belt objects, are remnants of the planetesimals that formed through this process during the formation of the Solar System about 4.5 billion years ago. Therefore, the minor body populations provide a glimpse into the past and offer a unique possibility to study the properties of planetesimals. Minor bodies, especially the ones that are not massive enough for self-gravity to press them into spheres, often have irregular shapes which might be reminiscent of their formation or caused by evolutionary processes throughout the history of the Solar System.

The main asteroid belt between the orbits of Mars and Jupiter is a collisionally evolved minor body population and only objects with diameters larger than ~ 120 km are thought to be primordial, which means that the properties of the large objects were shaped during the accretion process, while the smaller asteroids are by-products of collisional fragmentation (Bottke et al. 2005). The largest asteroids, for example Ceres and Vesta, are massive enough to be spherical. Smaller asteroids are rubble piles that formed through gravitational re-accumulation of the fragments from catastrophic collisions of larger bodies (Michel & Richardson 2013; Michel et al. 2020). Therefore, the irregular shapes of kilometre-sized asteroids are not primordial and not related to the planetesimal formation process.

Comets formed outside the water ice line where, besides refractory dust, volatiles were also present as ices. The favoured formation scenario for comets is the gravitational collapse of a pebble cloud, which is consistent with the observed properties of cometary nuclei (Blum et al. 2014, 2017, 2022). From spacecraft flybys it is known that cometary nuclei are irregularly shaped, many are elongated or bilobed, with sizes in the range of kilometres to a few tens of kilometres (Keller et al. 1986; Buratti et al. 2004; Duxbury et al. 2004; A’Hearn et al. 2005, 2011; Sierks et al. 2015). When entering the inner Solar System, solar irradiation triggers the sublimation of volatile ices and activity-driven erosion changes the surface morphology of the nuclei. Anisotropic mass loss caused by non-uniform insolation of the nucleus is capable of producing irregular-shaped nuclei (Vavilov et al. 2019). While transitioning from the comet reservoir to the inner Solar System, sublimation of volatiles can spin up the nucleus to disruption and the subsequent reassembly of the fragments can also lead to a bilobed nucleus (Safrit et al. 2021). Furthermore, modelling the collisional history of comet-sized planetesimals during the dynamical evolution of the outer Solar System suggests that kilometre-sized cometary nuclei are not primordial but predominantly the fragments of larger bodies (Morbidelli & Rickman 2015), although the collision history depends very strongly on the unknown lifetime of the primordial, massive Kuiper belt before the outward migration of Neptune. Numerical simulations have demonstrated that the re-accumulation of the fragments produced in sub-catastrophic or catastrophic collisions would lead to the elongated or bilobed shapes of cometary nuclei while preserving volatiles and porosity (Jutzi & Benz 2017; Schwartz et al. 2018). This shows that the shapes of cometary nuclei as observed today are not necessarily primordial, but can originate from evolutionary processes as well.

Kuiper belt objects (KBOs) orbiting in the region beyond Neptune are another population of minor bodies in the Solar System. The inclination distribution of the Kuiper belt objects shows that there are two populations, a hot population with a wide inclination distribution and a cold population with a narrow inclination distribution (Brown 2001). The transition between both populations is set around 5° and KBOs that have an inclination ≲5° belong more likely to the cold population. This cold population, the so-called cold classical Kuiper belt, forms a narrow ring between approximately 42.4 au and 47.7 au, with a sub-population centred at ~44 au (the so-called kernel), and consists most likely of planetesimals that formed in situ, unaffected by the dynamical history of the outer Solar System formation (Batygin et al. 2011; Nesvorný et al. 2011; Petit et al. 2011; Kavelaars et al. 2021).

Kavelaars et al. (2021) compared the size distributions of the cold classical Kuiper belt with the size distributions typically found in streaming instability simulations. They find a good match with the exponentially tapered power law proposed in Schäfer et al. (2017). However, the total mass of the classical Kuiper belt is at least a factor 10 lower than the result of streaming instability; this may be due to gravitational perturbations from Neptune (Gomes et al. 2018). We refer readers to Simon et al. (2022) for a detailed discussion of the streaming instability initial-mass function compared to the classical Kuiper belt.

Observations show that the binary fraction among the cold classical Kuiper belt objects is very high (Noll et al. 2008a,b) and it is even hypothesised that almost all objects of the Kuiper belt formed as binaries (Fraser et al. 2017). The high binary fraction and the mostly prograde mutual orbits of the binaries are in agreement with high-resolution studies of the gravitational collapse of pebble clouds (Nesvorný et al. 2010, 2019, 2021; Robinson et al. 2020; Polak & Klahr 2023). Light curve modelling of KBOs revealed that besides wide binaries, contact binaries are also present among the KBOs. From the first confirmed contact binary 2001 QG298, a contact-binary fraction of ≳ 10–30% was estimated (Sheppard & Jewitt 2004; Lacerda 2011). Other candidates are the Kuiper belt objects 2003 SQ317 and 2004 TT357 that have light curves consistent with those of contact binaries even though highly elongated shapes cannot be entirely ruled out (Lacerda et al. 2014; Thirouin et al. 2017). The contact-binary fraction among the Plutinos, which are KBOs such as Pluto that are in a 3:2 resonance with Neptune, is even higher with an estimate of ~50% based on observations (Thirouin & Sheppard 2018). This shows that contact binaries are frequent with a fraction of 10–25% among the cold classical Kuiper belt objects and ~50% among the Plutinos (Thirouin & Sheppard 2019; Showalter et al. 2021). The cold classical Kuiper belt object (486958) Arrokoth is a contact binary that was visited by the National Aeronautics and Space Administration (NASA) spacecraft New Horizons (Stern et al. 2019). As a member of the cold classical Kuiper belt and a contact binary, Arrokoth most likely represents a primordial planetesimal that has formed by the gravitational collapse of a pebble cloud (Stern et al. 2019; McKinnon et al. 2020). Even though it is hypothesised that the shape of Arrokoth might be the result of volatile sublimation on a 1–100 Myr timescale (Zhao et al. 2021), its shape might still provide useful insights in the primordial shapes of planetesimals.

While the small asteroids are the results of catastrophic collisions and gravitational re-accumulation, comets and the cold classical Kuiper belt objects might still be representative of the shapes of the kilometre-sized primordial planetesimals that formed through the gravitational collapse of pebble clouds. The dynamics of the pebble cloud collapse has been studied in high resolution N-body simulations or by means of hydrodynamic simulations of the self-gravitating pebble fluid (Nesvorný et al. 2010, 2019, 2021; Robinson et al. 2020; Polak & Klahr 2023). However, these studies do not resolve the final shapes of the planetesimals. Here, we attempt to model the final assembly of a kilometre-sized planetesimal from the gravitational collapse of a rotating clump of pebbles using a Monte-Carlo method for granular dynamics and show that, depending on the angular momentum content, the final shapes of kilometre-sized planetesimals are in agreement with the shapes of the left-over planetesimals in the Solar System.

The paper is organised as follows. In Sect. 2, we broadly summarise available shape information of minor bodies. In Sect. 3, we provide an overview of the Monte-Carlo method, with a more detailed description in the Appendix A. In Sect. 4, we go throught the results of our simulations. In Sect. 5, we discuss our results in the context of the Solar System bodies. Finally, in Sect. 6, we summarise the main findings of our study.

2 Shapes of minor bodies

In this section, we present a brief overview of minor bodies with well-studied sizes, namely the small number of comets that were visited by space craft, the cold classical Kuiper belt object Arrokoth, and the interstellar asteroid 11/2017 U1 (‘Oumuamua).

2.1 Arrokoth

The overall dimensions of Arrokoth are 35.95 × 19.90 × 9.75 km. Because the flyby of the New Horizons spacecraft was fast and distant and because Arrokoth lacks any natural satellites, there are no direct constraints on the mass or the density by gravity measurements (Keane et al. 2022). The density is hence inferred from the geophysical environment of Arrokoth which favours a nominal density of ~235 kg m−3 and yields a mass of −7.485 × 1014 kg for the entire body (Keane et al. 2022). The two planetesimals that form the lobes of Arrokoth have dimensions of 21.20 × 19.90 × 9.05 km for the large lobe Wenu and 15.75 × 13.85 × 9.75 km for the small lobe Weeyo (Keane et al. 2022). The large lobe is hence significantly flattened along one axis, while the small lobe is slightly less flat. It is hypothesised that the flattening could be the result of volatile sublimation on a 1–100 Myr timescale (Zhao et al. 2021). With a bulk density of ~235 kg m−3, Arrokoth is less dense than other minor bodies of the Solar System and especially less dense than comets that have an estimated mean density of around 500 kg m−3. For a typical grain density of 2154 kg m−3 (rock mass-fraction of 75% with an average rock density of 3500 kg m−3), the volume-filling fraction of Arrokoth is only 11%, which translates to a porosity of 89% (Keane et al. 2022). As a member of the cold classical Kuiper belt and a contact binary, Arrokoth most likely represents a primordial planetesimal that has formed through the gravitational collapse of a pebble cloud and the subsequent low-velocity merger of two components (Stern et al. 2019; McKinnon et al. 2020; Marohnic et al. 2021). Therefore, the shape of Arrokoth provides useful insights in the primordial shapes of planetesimals.

2.2 Comets

Because of their small sizes, cometary nuclei are difficult to observe from Earth. However, information about their sizes and shapes can be obtained through light curve analysis or radar observation. Additionally, for six comets that were visited by spacecraft, detailed information about the shapes of their nuclei exists. The Giotto mission of the European Space Agency (ESA) to comet 1P/Halley obtained for the first time a direct view of a cometary nucleus and revealed an irregular-shaped kilometre-sized object (Keller et al. 1986, 1987). Further NASA missions to comets 19P/Borrelly (Deep Space 1), 9P/Tempel 1 and 103P/Hartley 2 (Deep Impact), and 81P/Wild 2 (Stardust) provided more examples for the irregular shapes of cometary nuclei (Buratti et al. 2004; Duxbury et al. 2004; A’Hearn et al. 2005, 2011). Lastly, the most recent mission to a comet, ESA’s Rosetta mission, which accompanied the Jupiter-family comet 67P/Churyumov-Gerasimenko for almost 2 yr along its orbit, obtained high-resolution images and detailed shape information of a kilometre-sized bilobed nucleus that is thought to be a contact binary of two distinct lobes (Massironi et al. 2015; Nesvorný et al. 2018), with a large lobe of size 4.1 × 3.5 × 1.6 km and a small lobe of size 2.5 × 2.1 × 1.6 km (Sierks et al. 2015; Jorda et al. 2016). The two lobes are not spherical but flattened along one axis. The nucleus of 67P has a mass of 9.982 × 1012 kg, a density of 532 kg m−3, and a high porosity of the order of 70% (Jorda et al. 2016). Shape models are also available for comets 103P/Hartley 2 and 19P/Borrelly. Both comets are bilobed and it is speculated that this is the result of formation from individual objects (Oberst et al. 2004; Thomas et al. 2013). Radar observations of comet 8P/Tuttle revealed a contact binary with a large lobe of size 5.75 × 4.11 km and a small lobe of size 4.25 × 3.27 km (Harmon et al. 2010).

2.3 1l/2017 U1 (‘Oumuamua)

The object 1I/2017 U1 (‘Oumuamua) was the first detected interstellar object to pass through the Solar System (Meech et al. 2017). Light curve analysis and shape modelling revealed that ‘Oumuamua is unusually elongated with an axes ratio as high as 6:1 suggesting an ellipsoid body with axes diameters of 230 × 35 × 35 m (Jewitt et al. 2017). The overall best-fit model of Mashchenko (2019) suggests a thin disc with dimensions 115 × 111 × 19 m. A cigar-shaped body with dimensions of 324 × 42 × 42 m is also possible, but less likely according to the analysis of Mashchenko (2019), as the thin disc model fits the observed light-curve minima best.

3 Methods

To simulate the shape of a planetesimals, it is necessary to study the gravitational collapse of a pebble cloud. However, even a kilometre-sized planetesimal requires ~1018 millimetre-sized pebbles to form. It is hence computationally impossible to do this directly by following each pebble individually. However, one can simulate the time evolution of the distribution function of the pebbles subject to collisions and self-gravity, or in other words, solve the Boltzmann equation. We do this here by using a direct-simulation Monte Carlo (DSMC) method (Bird 1994; Alexander & Garcia 1997; Pöschel & Schwager 2005). The key idea of DSMC is to sample the distribution function with a number of computational particles, each of which represents a swarm of physical pebbles with equal properties. The Boltzmann equation, which governs the time evolution of the distribution function, is then split into an advection step, in which the positions of the particles are propagated with their current velocities, and a collision step, in which the particle velocities change due to collisions (see Appendix A for a detailed description). Self-gravity is accounted for in the advection step by using a drift-kick-drift scheme and a Poisson solver based on the Fast Fourier Transformation (FFT) (see Appendix B).

The collision treatment of the DSMC method based on the Boltzmann equation was originally developed to simulate dilute gases. In this case, the assumption is that two particles move independently from each other and that the two-particle distribution function can be written as the product of two one-particle distribution functions f2(x1,x2,v1,v2,t)=f(x1,v1,t)f(x2,v2,t)${f_2}\left( {{{\bf{x}}_1},{{\bf{x}}_2},{{\bf{v}}_1},{{\bf{v}}_2},t} \right) = f\left( {{{\bf{x}}_1},{{\bf{v}}_1},t} \right)f\left( {{{\bf{x}}_2},{{\bf{v}}_2},t} \right)$(1)

(molecular chaos hypothesis). For dense gases, however, this assumption may not hold any more because of correlations between the particles due to their finite volumes. The most simple correction to account for the finite volume effects is to introduce the Enskog factor χ, which is basically the two-particle correlation function, such that the distribution function reads f2(x1,x2,v1,v2,t)=χf(x1,v1,t)f(x2,v2,t)${f_2}\left( {{{\bf{x}}_1},{{\bf{x}}_2},{{\bf{v}}_1},{{\bf{v}}_2},t} \right) = \chi f\left( {{{\bf{x}}_1},{{\bf{v}}_1},t} \right)f\left( {{{\bf{x}}_2},{{\bf{v}}_2},t} \right).$(2)

(Pöschel & Schwager 2005). The Enskog factor χ depends on the local density of particles and can be derived from the equation of state of a hard-sphere fluid (Ma & Ahmadi 1986; see Appendix A). The Boltzmann equation then turns into the Boltzmann-Enskog equation, which is the basis for our study (Montanero & Santos 1996, 1997).

To spatially resolve the planetesimal, space is discretised using a Cartesian grid. While the Enskog correction makes sure that the collision rates of the DSMC method are correct, it does not affect the propagation of particles. In spatially resolved simulations, particles hence move unconditionally and unaffected by other particles, which can lead to an unrealistically high packing of particles within cells. To prevent this, the propagation step is modified and new positions are accepted based on the local packing of particles (Hong & Morris 2021). This probabilistic approach prevents the over-packing of cells and further allows to identify fully packed cells as the solid elements comprising the planetesimal. The maximum volume-filling factor of a cell is ~0.64, which is achieved for a random close packing of hard spheres.

A significant limitation of the method is that because of the fixed grid, it is not possible to properly resolve rotation and hence binary formation as once cells fill up, the bulk of the particles remains in the cell as the mean-free path is so short that only particles close to the cell boundaries might be able to diffuse outwards. However, for clumps of low angular momentum that collapse into a single body, the method is able to reproduce the shapes of planetesimals.

4 Simulating the final assembly

We focus on the final assembly of a kilometre-sized planetesimal from a pebble cloud with low angular momentum content. We cannot resolve the entire collapse from a Hill-radius-sized pebble cloud down to a planetesimal because the required number of grid cells to achieve sub-kilometre resolution would be too high. Instead, we assume that the pebble cloud has already fragmented into sub-clumps that would eventually form, for example, the components of a binary as has been shown in N-body simulations of the gravitational collapse (e.g. Polak & Klahr 2023).

thumbnail Fig. 1

Shapes of planetesimals for different values of the initial angular momentum. The rows show the column densities of cells that have a volume-filling factor ≥0.5 in the xy- and xz-planes for different values of the initial angular momentum of the cloud at the end of each simulation. We omit the yz-plane because of the evident cylindrical symmetry. The contour of the ellipsoid fitted to the planetesimals in the respective plane are shown in black with the centre of the ellipse marked with a black circle. The figure shows the results for initial values of L/LJ = 1.5, 1.0, 0.5 and 0.0 to illustrate the effect of angular momentum of the cloud on the flattening of the planetesimals.

4.1 Initial conditions

We started all simulations with an initial radius of the pebble cloud that was 10 times larger than the spherical planetesimal size corresponding to the total mass. The simulation box and the number of grid cells were chosen such that we achieved a resolution of 1/10 of the planetesimal radius. For a planetesimal of radius Rp = 1 km (as chosen in this study) this meant that the initial cloud radius was Rc = 10 km and the resolution was 0.1 km. The pebbles had a fixed radius of 0.5 mm in agreement with expected pebble sizes (Zsom et al. 2010; Lorek et al. 2018) and a constant coefficient of restitution of 0.5 (Weidling & Blum 2015). The pebbles were initially uniformly distributed within the cloud volume and their initial velocities were isotropic following a Maxwellian distribution with kBT/m≈(1/3)υvir$\sqrt {{k_{\rm{B}}}T/m} \approx (1/3){\upsilon _{{\rm{vir}}}}$, where υvir is the virial velocity of the pebble cloud. In addition, we added rotation to the pebble cloud by setting an initial angular momentum L. A critically rotating Jacobi ellipsoid of mass Mp and effective radius Rp has an angular momentum of LJ=0.4GMp3Rp.${L_{\rm{J}}} = 0.4\sqrt {GM_{\rm{p}}^3{R_{\rm{p}}}} .$(3)

We therefore set the angular momentum content in units of LJ because the cloud can collapse into a single body only for L/LJ ≲ 1. To study the effect of angular momentum on the planetesimal shape we systematically varied L/LJ between 0 (no rotation) and 2 (super-critical rotation).

We conducted the study for a pebble cloud that had the mass of a planetesimal of Rp = 1 km and density 1280kg m−3, which corresponded to a volume-filling factor of the planetesimal of 0.64 and a material density of the pebbles of 2000 kg m3. The pebble cloud was considered in isolation and therefore its dynamics was independent of the heliocentric distance and rotation around the Sun. This is a reasonable assumption because the free-fall time for such a small cloud is of the order of days, which is much shorter than the orbital period in the Kuiper belt, and furthermore, the size of the cloud is significantly smaller than the Hill radius so that shearing effects can be neglected.

To identify a planetesimal, we ran the simulations for typically ~ 1.2–1.4 free-fall times of the cloud until a solid body formed that no longer changed its shape significantly. The free-fall time was of the order of a day. We then identified all cells that had a volume-filling factor ≥0.5, which is the loosest stable packing of spheres (Torquato & Stillinger 2007). Because we considered the pebbles as spheres and because we did not have a size distribution of pebbles, we set the maximum volume-filling factor that a cell could obtain to 0.64, which is obtained for a random close packing of spheres and slightly higher than that for a random loose packing of spheres in the zero-gravity limit with volume-filling factor of 0.56 (Onoda & Liniger 1990).

4.2 Planetesimal shapes

Figure 1 shows a selection of planetesimals that formed from clouds with different initial angular momentum. With increasing angular momentum of the cloud, the planetesimals flatten in the z-direction which is parallel to the initial rotation axis. In the plane perpendicular to the rotation axis (xy-plane), the planetesimals remain more round. As expected for very high angular momentum content L/LJ ≳ 1.5, the cloud does not collapse into a single object but instead forms a flat disc with spiral structures.

Those clouds would possibly form binary systems. For lower angular momentum in the range 1 ≲ L/LJ ≲ 1.5, very thin disclike objects form. When the angular momentum is L/LJ ≲ 1, the planetesimal obtain lenticular and nearly spherical shapes comparable to the shapes of ‘Oumuamua or the lobes of Arrokoth and 67P (see discussion in Sect. 5). In summary, we find disc-like and lenticular-shaped planetesimals in our simulations. However, we do not see elongated ellipsoids that would resemble bodies such as 19P/Borrelly or 103P/Hartley 2.

Table 1

Comet and KBO sizes.

4.3 Axes ratios of planetesimals

We used linear least-squares to fit an ellipsoid that is centred at the centre-of-mass of the planetesimal to the boundary points (see Appendix C). The orientation of the ellipsoid does not necessarily coincide with the Cartesian axes and we hence determined the eigenvalues and eigenvectors of the ellipsoid which give the lengths and the directions of the main axes. Because of the limited resolution of the Cartesian grid, the ellipsoid is only an approximation to quantify the (irregular) shape of the planetesimal. To check the validity of the approximation, however, we calculated the inertia tensor of the planetesimal and compared the directions of the principal axes with the main axes of the ellipsoid. The axes were aligned, which showed that the ellipsoid orientation matched the mass distribution of the body. We quantify the shape of the planetesimals by calculating the ratios of the three axes of the fitted ellipsoid. We labelled the axes a, b, and c such that c ≤ b ≤ a and calculated the ratios b/a, c/a, and c/b. The closer any of the ratios is to unity, the more spherical is the object. For example, a planetesimal with b/a = 1, but c/a = 0.5 and c/b = 0.5, would be round in the ba-plane and flattened in the ca- and cb-planes resulting in a lenticular shape.

Figure 2 shows the axes ratios of the planetesimals as a function of the initial angular momentum content. The figure confirms that for L/LJ ≲ 1 the planetesimals become more spherical as is already indicated in Fig. 1. For comparison, we plot the axes ratios for the two lobes of Arrokoth (see Table 1 for the shape information). Our model is consistent with the axes ratios, and hence the shapes of the two planetesimals forming Arrokoth, for pebble clouds with angular momentum approximately in the range 0.5 ≲ L/LJ ≲ 0.9. Here, we chose Arrokoth for comparison because being a cold classical Kuiper belt object, Arrokoth most likely represents a primordial planetesimal. Figure 3 shows the average axes ratios for planetesimals that formed from clumps with initial angular momentum L/LJ in the range 0.5–0.9. For each value of L/LJ, we averaged the axes ratios over five realisations for the same initial conditions to obtain a measure for the uncertainty of our Monte-Carlo simulations. Figure 3 shows that the uncertainty of the axes ratios are ~10%.

thumbnail Fig. 2

Axes ratios of planetesimals as a function of the initial angular momentum of the cloud. The axes ratios of the two lobes of Arrokoth are shown as dashed (big lobe) and dotted (small lobe) lines.
thumbnail Fig. 3

Average axes ratios of planetesimals for the initial angular momentum of the cloud in the range 0.5 to 0.9. For each initial L/LJ. the results of 5 different realisations are averaged. The axes ratios of the two lobes of Arrokoth are shown as dashed (big lobe) and dotted (small lobe) lines.
thumbnail Fig. 4

Angular momentum of planetesimals as a function of the initial angular momentum of the cloud. The dashed line indicates where the final angular momentum of the planetesimal Lfinal equals the initial angular momentum and the dotted line where Lfinal equals half of the initial value. The colour coding corresponds to the mass that ends up in the planetesimal. The red crosses show the total angular momentum Llot in the simulation.

4.4 Angular momentum

Figure 4 shows the final angular momentum of the planetesimals as a function of the initial angular momentum of the pebble cloud. We find that the final planetesimal contains about 50% of the initial angular momentum and that between 60–80% of the mass of the pebble cloud ends up in the planetesimal. The figure also shows that the total angular momentum is conserved within ~5%. This deviation is the result of the open boundary conditions and the loss of particles that leave the simulation box taking away a tiny fraction of the total angular momentum.

5 Discussion

In this section, we discuss our results in the context of the minor bodies of the Solar System. We provided a brief overview of the processes that shape various minor bodies in the Solar System in Sect. 1. Small asteroids are typically rubble piles that formed through the gravitational reassembly of fragments produced in catastrophic collisions of larger bodies (Michel & Richardson 2013; Michel et al. 2020). The shapes of kilometre-sized asteroids are hence most likely not primordial and not related to the gravitational collapse of pebble clouds. We therefore exclude asteroids from our discussion. Cometary nuclei have irregular shapes that can originate in their formation process or be the result of evolutionary process (e.g. Jutzi & Asphaug 2015; Jutzi & Benz 2017; Schwartz et al. 2018; Safrit et al. 2021). Especially the shapes of the comets that are clearly contact binaries might provide some information about the formation in collapsing pebble clouds. Lastly, we include Arrokoth in our discussion. As a cold classical Kuiper belt objects, Arrokoth resembles best a primordial planetesimal (Stern et al. 2019; McKinnon et al. 2020). Table 1 summarises the available shape information for Solar System bodies that we use in our comparison.

5.1 Comparison with simulations

Figure 5 visualises the axes ratios of Solar System bodies, ‘Oumuamua, and the planetesimals formed in our simulations. It is clear from this comparison that most of the Solar System bodies have axes ratios different to what we find in the simulations. The planetesimals formed in our simulations are very close to oblate spheroids with b/a ≈ 1 and c/a ≲ 1. The higher the angular momentum of the pebble cloud, the more flattened are the planetesimals. The disc-like shape of ‘Oumuamua would be consistent with the very flat planetesimals that form in pebble clouds with 1 ≲ L/LJ ≲ 1.5. On the other hand, most comets are prolate in shape or triaxial ellipsoids with c ≈ b ≲ a (see also Table 1). Arrokoth as a whole has different axes ratios. Therefore, the observed shapes do not seem to be in agreement with the shapes of planetesimals formed by gravitational collapse. However, we neglected the fact that four out of the six comets that were visited by spacecraft are bilobed (19P/Borrelly, 67P/C-G, 103P/Hartley 2) or potentially bilobed (1P/Halley). Among the bilobed comets, 67P is a contact binary with two distinct lobes for which shape information is available (Sierks et al. 2015; Jorda et al. 2016). Furthermore, comet 8P/Tuttle is a contact binary as confirmed by radar observations (Harmon et al. 2010). And lastly, the cold classical Kuiper belt object Arrokoth is a contact binary as well with well characterised shape (Keane et al. 2022). Comets 19P/Borrelly and 103P/Hartley 2 are bilobed objects that might be the result of the formation from two individual objects (Britt et al. 2004; Oberst et al. 2004; Thomas et al. 2013). Using the available shape model for 103P/Hartley 2 (Farnham & Thomas 2013), we fitted ellipsoids to the two lobes of the comet (see Appendix D). Comet 19P/Borrelly lacks a detailed shape model and we could not extract the two lobes. For 1Ρ/Halley it is not entirely clear whether the nucleus is bilobed or not.

As outlined in Sect. 1, planetesimals are thought to form by gravitational collapse of pebble clouds. The observed properties of comets and the characteristics of Kuiper belt objects, for example the size distribution and the occurrence of binaries, are consistent with this formation hypothesis (Blum et al. 2014, 2017, 2022; Nesvorný et al. 2010, 2019, 2021; Polak & Klahr 2023). The dynamics of the collapse has been studied in detail with N-body simulations and shows that the pebble cloud typically fragments into a number of sub-clumps that eventually form planetesimals and bound binary systems (Nesvorný et al. 2010, 2021; Robinson et al. 2020; Polak & Klahr 2023). The subsequent collapse of the binary or low-velocity collisions of two planetesimals would lead to the formation of contact binaries, such as Arrokoth or 67P. In this framework, the lobes of the contact binaries would be planetesimals that formed through gravitational collapse of the sub-clumps. Therefore, it is more reasonable to compare the shapes of the planetesimals in our simulations with the shapes of the lobes of the contact binaries. The right panel of Fig. 5 shows how the simulated planetesimals compare to the lobes of the contact binaries. Besides the large lobes of 8P and 103P, the axes ratios of the other bodies overlap with the simulated planetesimals. This strongly supports the notion that the planetesimals forming a contact binary form from the collapse of a pebble cloud where angular momentum and gravity create oblate spheroids.

thumbnail Fig. 5

Axes ratios of Solar System bodies, ‘Oumuamua, and simulated planetesimals. The simulated planetesimals are shown as circles with a colour that indicates the initial angular momentum of the pebble cloud L/LJ. For L/LJ > 1, we used open symbols. The data points with error bars are from simulations where we averaged over five runs with varying random seed. The Solar System bodies and ‘Oumuamua are shown with open diamond symbols and different colours: 1P/Halley (blue), 8P/Tuttle (orange), 9P/Tempel 1 (green), 19P/Borrelly (red), 67P/Churyumov-Gerasimenko (purple), 81P/Wild 2 (brown), 103P/Hartley 2 (pink), Arrokoth (grey), 1I/2017 U1 (‘Oumuamua) (olive). Left: axes ratios of the whole bodies. Right: axes ratios of the individual lobes of the bilobed and contact-binary comets, Arrokoth, and the simulated planetesimals. We plot the axes ratios of the individual lobes with an upper triangle for the big lobe and a lower triangle for the small lobe. For Hartley 2, we included the two lobes that we found from analysing the shape model (see Appendix D).

5.2 The uncertain origin of 1I/2017 U1 (‘Oumuamua)

While the formation and nature of ‘Oumuamua is still under debate, a natural origin is advocated (‘Oumuamua ISSI Team et al. 2019). Different formation hypotheses are discussed in the literature that build on the observations of ‘Oumuamua, such as the unusually elongated shape (Jewitt et al. 2017; Meech et al. 2017), the lack of visible cometary activity (Jewitt et al. 2017; Meech et al. 2017; Trilling et al. 2018), and the non-gravitational acceleration (Micheli et al. 2018). It is unlikely that the non-gravitational acceleration is caused by strong volatile out-gassing, as it is the case for Solar System comets, because the outgassing torques would spin up ‘Oumuamua to rotational fission (Rafikov 2018).

Flekkøy et al. (2019) discuss the possibility of ‘Oumuamua being a fractal dust aggregate, which would be consistent with their estimates of the mechanical stability and the change in rotation period due to radiation pressure. Seligman & Laughlin (2020) hypothesise that ‘Oumuamua could be a body rich in molecular hydrogen ice that formed in a cold dense core of a giant molecular cloud and that the sublimation of H2 could explain the non-gravitational acceleration. Desch & Jackson (2021, 2022) and Jackson & Desch (2021) argue that ‘Oumuamua could be a nitrogen-ice fragment excavated from an exo-Pluto with its shape heavily affected by N2-ice sublimation while approaching the Sun. Raymond et al. (2018) use dynamical simulations to show that ‘Oumuamua could be an extinct fragment of an ejected cometary planetesimal. Bergner & Seligman (2023) interpret ‘Oumuamua as an icy comet-like planetesimal. They show that the non-gravitational acceleration of ‘Oumuamua could be driven by the outgassing of molecular hydrogen that was created by cosmic rays, trapped in the water ice matrix, and released upon annealing of the amorphous water ice matrix while approaching the Sun. Farnocchia et al. (2023) and Seligman et al. (2023) report the detection of sub-kilometre-sized near-Earth asteroids that exhibit significant non-gravitational acceleration that is greater than what can be attributed to the Yarkovsky effect. Similar to ‘Omuamua, these objects do not show any sign of visible coma, but weak volatile outgassing could explain their orbital motion.

In summary, despite the different hypotheses, the origin of ‘Oumuamua remains elusive. Being agnostic about the actual origin of 'Oumuamua, our results demonstrate that the disc-like flattened shape of 'Oumuamua could arise naturally from the gravitational collapse of a pebble cloud with high angular momentum.

5.3 Limitations of the method

Our method is capable of resolving the shapes of planetesimals formed by gravitational collapse. However, the fixed grid that is necessary for resolving pebble collisions sets some limitations for our model. In order to resolve the substructure of kilometre-sized planetesimal, we cannot follow the entire collapse from a Hill-sized pebble cloud to a planetesimal as it would require too many grid cells. Therefore, we are limited to follow the final assembly of a planetesimal from sub-clumps that form during the collapse (Nesvorný et al. 2021; Polak & Klahr 2023). Because of the fixed grid, cells that are filled to the maximum volume-filling factor will not move any more. Therefore, we cannot properly resolve the rotational state of the final planetesimal. A similar problem arises for pebble clouds with high angular momentum that collapse into disc-like flat structures (L/LJ ≳ 1). Once the disc structure has formed, the dense cells are frozen out. In reality, gravity and angular momentum would cause the substructures to collapse and form binaries or contact binaries. Therefore, our method is reliable for low angular-momentum clouds that collapse into a single body.

6 Conclusion

We developed a Monte-Carlo method to simulate the collapse of a rotating pebble cloud and the final assembly of a planetesimal. We investigated the shapes by analysing the axes ratios of ellipsoids fitted to the final body and compared our results to the shapes of Solar System comets with known shapes, the interstellar asteroid 'Oumuamua, and the cold classical Kuiper belt object Arrokoth. We find that planetesimals that form by gravitational collapse of pebble clouds contain about 60–80% of the mass of the cloud and ~50% of the initial angular momentum. Furthermore, planetesimals that form from pebble clouds with initial angular momentum in the range 0.5 ≲ L/LJ ≲ 1, where LJ is the angular momentum of a critically rotating Jacobi ellipsoid with the same mass and effective radius as the planetesimal, match the observed shapes of the lobes of Arrokoth. The shapes of the lobes of bilobed comets agree with the shapes of the planetesimals that form from L/LJ ≲ 1. Lastly, the disc-like shape of the interstellar asteroid 'Oumuamua, if this body is indeed primordial and not a collisional fragment, is consistent with a formation in pebble clouds with 1 ≲ L/LJ ≲ 1.5. Because comets and, especially, Arrokoth most likely represent the primordial planetesimals of the Solar System (even though evolutionary processes might have lead to some reshaping), our results suggest that the flattened shapes of planetesimals naturally follows from the gravitational collapse of a rotating pebble cloud. This further strengthens the hypothesis of planetesimal formation by gravitational collapse.

Acknowledgements

We thank the anonymous referee for helpful comments that helped improving this manuscript. A.J. acknowledges funding from the European Research Foundation (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 2018-04867), the Danish National Research Foundation (DNRF Chair Grant DNRF159) and the Goran Gustafsson Foundation.

Appendix A Monte Carlo method for solving the Enskog equation

The Monte Carlo method for solving the Enskog equation follows the procedure outlined in Montanero & Santos (1997) and Pöschel & Schwager (2005). We use a total of N particles to sample the distribution function of the pebbles. Each particle represents a swarm of Nswm = Nphys/N pebbles, where Nphys =Mp/mpeb is the total number of pebbles of mass mpeb that is needed for a planetesimal of mass Mp. Each particle has position x, velocity v, mass m, and radius R. The Monte-Carlo method calculates the change of particle velocities owing to collisions.

The number of collisions of particle i in cell I with particles j in cell J in a time step Δt is ωij=σ24π(vij⋅ei)χ(xi,xi+σei)nJΔt.${\omega _{ij}} = {\sigma ^2}4\pi \left( {{{\bf{v}}_{ij}} \cdot {{\bf{e}}_i}} \right)\chi \left( {{{\bf{x}}_i},{{\bf{x}}_i} + \sigma {{\bf{e}}_i}} \right){n_J}\Delta t.$(A.1)

Here, σ is the diameter of a pebble in case of a mono-disperse size distribution of pebbles; for a non-uniform size distribution, σ is chosen at random in the interval between the minimum and the maximum pebble radius such that Ri + Rmin≤σ≤Ri+Rmax. The relative velocity between particles i and j is vij = vj − vi and ei is the direction of the collision, that is, the unit vector pointing from particle i to j. Finally, χ(xi, xi + σei) is the local equilibrium pair correlation function (or Enskog factor) (Van Beijeren & Ernst 1973), which depends on the density in the cell J, and nJ is the number density in cell J. The local equilibrium pair correlation function can be determined from the equation of state (Frezzotti 1998). Here, we use the hard-sphere equation of state from Ma & Ahmadi (1986), which gives χ=1+2.5ϕ+4.5904ϕ2+4.515439ϕ3(1−(ϕ/ϕmax)3)0.67802$\chi = {{1 + 2.5\phi + 4.5904{\phi ^2} + 4.515439{\phi ^3}} \over {{{\left( {1 - {{\left( {\phi /{\phi _{\max }}} \right)}^3}} \right)}^{0.67802}}}}$(A.2)

for the correlation function, where ϕ is the packing fraction and ϕmax is the maximum packing fraction. The maximum packing fraction, we set to 0.6, which lies in the range between random loose packing (RLP) with ϕmax = 0.56 and random close packing (RCP) with ϕmax = 0.64 (Onoda & Liniger 1990).

The idea of the Monte Carlo method is to sample a certain number of collision for each cell. How many collisions are expected is calculated as Ncoll =12NIωmax, ${N_{{\rm{coll }}}} = {1 \over 2}{N_I}{\omega _{\max }}{\rm{, }}$(A.3)

where ωmax= maxij({ωij}), for any particle i of cell I and any particle j of cell J separated a distance σ from I, is an upper bound on the number of collisions (Montanero & Santos 1997). For each collision, we then carry out the following steps:

  • 1.

    A particle i from all particles in the current cell I is chosen at random with equiprobability.

  • 2.

    A direction ei is chosen at random with equiprobability.

  • 3.

    A particle j from the cell J where xi+σei points to is chosen at random.

  • 4.

    The collision between particles i and j is accepted with probability Θ(ei ⋅ vij)ωij/ωmax. Here, Θ(ei ⋅ vij) makes sure that the particle pair is approaching.

  • 5.

    If the collision is accepted, new velocities are assigned to particles i and j according to vi′=vi+(1+ϵ)mjmi+mjvij⋅ei and vj′=vj−(1+ϵ)mimi+mjvij⋅ei${\bf{v}}_i^\prime = {{\bf{v}}_i} + (1 + ){{{m_j}} \over {{m_i} + {m_j}}}{{\bf{v}}_{ij}} \cdot {{\bf{e}}_i}{\rm{ and }}{\bf{v}}_j^\prime = {{\bf{v}}_j} - (1 + ){{{m_i}} \over {{m_i} + {m_j}}}{{\bf{v}}_{ij}} \cdot {{\bf{e}}_i}$, where ϵ is the coefficient of restitution.

Once all collisions are sampled, we proceed with the next cell.

Even though the Enskog equation takes excluded volume and correlations between particles into account, it does not prevent cells from becoming unrealistically dense, that is, having ϕ > ϕmax. The reason for this is that when particles are propagated, their new positions are unconditionally accepted (Hong & Morris 2021). To prevent this from happening, new particle positions (if the particle enters a different cell) are accepted with a probability of exp[−6(χ(ϕ1)ϕ1 − χ(ϕ0)ϕ0] where ϕ0,1 and χ(ϕ0,1) are the packing fractions and correlation functions before and after the particle would enter the cell (Hong & Morris 2021).

Appendix B Self-gravity of the cloud

The Poisson equation, ΔΦ(x)=4πρ(x),$\Delta \Phi ({\bf{x}}) = 4\pi \rho ({\bf{x}})$(B.1)

determines the dynamical evolution of the pebble cloud under its own gravity. Here, Φ is the gravitational potential and ρ is the density of the pebble distribution. Poisson’s equation can be easily solved in Fourier space. The potential and the density are written in terms of their Fourier transforms Φ(x)=∫Φ^(k)ei2πk⋅xdk,$\Phi ({\bf{x}}) = \int {\hat \Phi } ({\bf{k}}){e^{i2\pi {\bf{k}} \cdot {\bf{x}}}}{\rm{d}}{\bf{k}},$(B.2) ρ(x)=∫ρ^(k)ei2πk⋅xdk,$\rho ({\bf{x}}) = \int {\hat \rho } ({\bf{k}}){e^{i2\pi {\bf{k}} \cdot {\bf{x}}}}{\rm{d}}{\bf{k}},$(B.3)

where k is the wave number and a hat denotes quantities in Fourier space. Poisson’s equation in Fourier space simplifies to the algebraic equation Φ^=−4πGρ^4π2|k|2$\hat \Phi = - 4\pi G{{\hat \rho } \over {4{\pi ^2}|{\bf{k}}{|^2}}}$(B.4)

We first calculate ρ^$\hat \rho $ using the Fast Fourier Transform (FFT), then solve the above equation, and in a last step calculate the inverse FFT of Φ^$\hat \Phi $ to obtain the potential Φ in real space. For the FFT, we make use of the FFTW package (Frigo & Johnson 2005), which is an highly optimised software package to carry out fast FFTs. We calculate the accelerations at the cell centres from the potential Φ using a standard central finite-difference scheme ax,ijk=−Φi+1jk−Φi−1jk2Δx,${a_{x,ijk}} = - {{{\Phi _{i + 1jk}} - {\Phi _{i - 1jk}}} \over {2\Delta x}},$(B.5) ay,ijk=−Φij+1k−Φij−1k2Δy, ${a_{y,ijk}} = - {{{\Phi _{ij + 1k}} - {\Phi _{ij - 1k}}} \over {2\Delta y}}{\rm{, }}$(B.6) az,ijk=−Φijk+1−Φijk−12Δz${a_{z,ijk}} = - {{{\Phi _{ijk + 1}} - {\Phi _{ijk - 1}}} \over {2\Delta z}}$(B.7)

where Δx, Δy, and Δz are the grid spacings in each direction. Finally, we perform a linear interpolation to obtain the acceleration at the actual particle position.

The Fourier method implicitly assumes periodic boundary conditions. Therefore, it is important to chose a simulation box large enough for boundary effects owing to the periodic images of the pebble cloud to be unimportant. We found from testing different box sizes for the pure gravitational collapse of the pebble cloud that a box which is about two times larger than the pebble cloud diameter is sufficient.

Appendix C Fitting an ellipsoid

To determine the axes ratios of the planetesimals, we fitted an ellipsoid to the data points. This can be done as we describe in the following. The general equation for a conic section is given by the polynomial Ax2+By2+Cz2+Dxy+Exz+Fyz+Gx+Hy+Iz+J=0.$A{x^2} + B{y^2} + C{z^2} + Dxy + Exz + Fyz + Gx + Hy + Iz + J = 0.$(C.1)

We have N data points (xi, yi, zi) for the boundary of the planetesimal, which provides us with a set of N equations of the above type of Eq. C.1. Normalising such that J = −1 and introducing the N × 9 matrix M which has rows formed by the vectors (xi2 yi2 zi2 xiyi xizi yizi xi yi zi)$\left( {x_i^2y_i^2z_i^2{x_i}{y_i}{x_i}{z_i}\,{y_i}{z_i}{x_i}{y_i}{z_i}} \right)$, we can write the matrix equation Mp=b, $M{\bf{p}} = {\bf{b}}{\rm{, }}$(C.2)

where we have further introduced two vectors, one of length 9 for the parameters p=(A B C D E F G H I)⊤ and one of length N for the right-hand side b=(1 … 1). The parameters of the ellipsoid are found by solving the above matrix equation C.2. Because M is not a square matrix, we multiply with M⊤ from the left side. The product M⊤ M is a square matrix which can be inverted to solve for p, which we obtain through p=(M⊤M)−1(M⊤b).${\bf{p}} = {\left( {{M^ \top }M} \right)^{ - 1}}\left( {{M^ \top }{\bf{b}}} \right).$(C.3)

Having the parameters, the next step is to extract the properties of the ellipsoid: shape, orientation, and centre. To do so, we first notice that the polynomial in Eq. C.1 can be written in matrix form as a quadratic form, y⊤Qy=0,${{\bf{y}}^ \top }Q{\bf{y}} = 0,$(C.4)

where we introduced the vector y=(x y z 1)⊤ and the matrix Q Q=(AD/2E/2G/2D/2BF/2H/2E/2F/2CI/2G/2H/2I/2J).$Q = \left( {\matrix{ A & {D/2} & {E/2} & {G/2} \cr {D/2} & B & {F/2} & {H/2} \cr {E/2} & {F/2} & C & {I/2} \cr {G/2} & {H/2} & {I/2} & J \cr } } \right).$(C.5)

To find the centre x0 = (x0 y0 z0) of the ellipsoid, we introduce a translation matrix of the form T=(100x0010y0001z00001)$T = \left( {\matrix{ 1 & 0 & 0 & {{x_0}} \cr 0 & 1 & 0 & {{y_0}} \cr 0 & 0 & 1 & {{z_0}} \cr 0 & 0 & 0 & 1 \cr } } \right)$(C.6)

such that Ty moves the ellipsoid to the centre x0. Doing so, we can write the shifted ellipsoid in matrix form as y⊤(T⊤ QT)y. The transformed matrix is symmetric and looks as follows T⊤QT=(AD/2E/2T14D/2BF/2T24E/2F/2CT34T14T24T34T44),${T^ \top }QT = \left( {\matrix{ A & {D/2} & {E/2} & {{T_{14}}} \cr {D/2} & B & {F/2} & {{T_{24}}} \cr {E/2} & {F/2} & C & {{T_{34}}} \cr {{T_{14}}} & {{T_{24}}} & {{T_{34}}} & {{T_{44}}} \cr } } \right),$(C.7)

where we notice that the upper left 3 × 3 matrix is the same as in Q. This is the part that describes the shape and orientation of the ellipsoid. The entries T14 to T34 arise from the fact that the ellipsoid is not centred at (0, 0, 0) but shifted to x0, that is, the terms linear in (x, y, z) in the polynomial form of the ellipsoid Eq. C.1 do not vanish. And finally, the entry T44 is the normalisation for recovering the enforced J = −1. The individual entries are T14=Ax0+D2y0+E2z0+G2,${T_{14}} = A{x_0} + {D \over 2}{y_0} + {E \over 2}{z_0} + {G \over 2},$(C.8) T24=D2x0+By0+F2z0+H2,${T_{24}} = {D \over 2}{x_0} + B{y_0} + {F \over 2}{z_0} + {H \over 2},$(C.9) T34=E2x0+F2y0+Cz0+I2,${T_{34}} = {E \over 2}{x_0} + {F \over 2}{y_0} + C{z_0} + {I \over 2},$(C.10) T44=Ax02+By02+Cz02+Dx0y0        +Ex0z0+Fy0z0+Gx0+Hy0+Iz0+J.$\eqalign{ & {T_{44}} = Ax_0^2 + By_0^2 + Cz_0^2 + D{x_0}{y_0} \cr & \,\,\,\,\,\,\,\, + E{x_0}{z_0} + F{y_0}{z_0} + G{x_0} + H{y_0} + I{z_0} + J. \cr} $(C.11)

To find the centre of the ellipsoid, we notice that the entries T14=T24=T34=0 need to be equal to zero. This gives the following linear matrix equation (AD/2E/2D/2BF/2E/2F/2C)(x0y0z0)=(−G/2−H/2−I/2), $\left( {\matrix{ A & {D/2} & {E/2} \cr {D/2} & B & {F/2} \cr {E/2} & {F/2} & C \cr } } \right)\left( {\matrix{ {{x_0}} \hfill \cr {{y_0}} \hfill \cr {{z_0}} \hfill \cr } } \right) = \left( {\matrix{ { - G/2} \hfill \cr { - H/2} \hfill \cr { - I/2} \hfill \cr } } \right){\rm{, }}$(C.12)

which can be easily solved to find the centre. Next we want to find the axes and orientation of the ellipsoid. We can do this by finding the eigenvalues and eigenvectors of the matrix ℰ=(AD/2E/2D/2BF/2E/2F/2C),${\cal E} = \left( {\matrix{ A & {D/2} & {E/2} \cr {D/2} & B & {F/2} \cr {E/2} & {F/2} & C \cr } } \right),$(C.13)

which we need to normalise by dividing all entries by −T44 in order to recover J = −1. The eigenvalues are the solution of the characteristic polynomial det(ℰ−λI3)=0,$\det \left( {{\cal E} - \lambda {_3}} \right) = 0,$(C.14)

where 𝕀3 is the 3 × 3 identity matrix and λ is the eigenvalue. Because in diagonal form the matrix describing the ellipsoid is ℰdiag =(λa000λb000λc),${{\cal E}_{{\rm{diag }}}} = \left( {\matrix{ {{\lambda _a}} & 0 & 0 \cr 0 & {{\lambda _b}} & 0 \cr 0 & 0 & {{\lambda _c}} \cr } } \right), $(C.15)

which results in λax2 + λby2 + λcz2 = 1, the axes are related to the eigenvalues as a=λa−1/2, b=λb−1/2$a = \lambda _a^{ - 1/2},\quad b = \lambda _b^{ - 1/2}$, and c=λc−1/2$c = \lambda _c^{ - 1/2}$. For each eigenvalue, a eigenvector q can be calculated by the standard procedure of solving the eigenvalue equation ℰq = λq, which gives the direction of the corresponding axis.

Appendix D Fitting two lobes to 103P/Hartley 2

Based on the shape model of 103P/Hartley 2 (Farnham & Thomas 2013), we cut the nucleus at the neck region in two parts. To do so, we identified the thinnest cross section of the neck and fitted a plane that separated the nucleus in a big and a small lobe, in a similar fashion as Keane et al. (2022) described it for Arrokoth. After that, we fitted ellipsoids to each lobe separately. The separation plane is given as −0.0088x+0.0866y+z=−0.4161, where (x, y, z) are planetocentric Cartesian coordinates. The fitting results are summarised in Table D.1 and we use our results for the big and small lobes in Fig. 5.

Table D.1

Ellipsoidal fit 103P/Hartley 2.

Appendix E 3D visualisation of a planetesimal

As an example, Fig. E.1 shows a 3D visualisation of a planetesimal formed in our simulations. The initial angular momentum of the cloud was L/LJ = 0.5. One can see the flattened structure of the planetesimal arising from angular momentum conservation.

thumbnail Fig. E.1

Shape model of a planetesimal with L/LJ=0.5. Each volume element with a volume-filling factor ≥0.5 is plotted to visualise the shape of the planetesimal. The centres of the volume elements are projected into the xy-, xz-, and yz-planes to visualise the shape of the planetesimal that can not be seen in the three-dimensional representation.

References

  1. A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258 [Google Scholar]
  2. A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396 [CrossRef] [Google Scholar]
  3. Alexander, F. J., & Garcia, A. L. 1997, Comput. Phys., 11, 588 [NASA ADS] [CrossRef] [Google Scholar]
  4. Batygin, K., Brown, M. E., & Fraser, W. C. 2011, ApJ, 738, 13 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bergner, J. B., & Seligman, D. Z. 2023, Nature, 615, 610 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bird, G. A. 1994, Molecular Gas Dynamics And The Direct Simulation Of Gas Flows (Oxford: Clarendon Press) [Google Scholar]
  7. Blum, J. 2018, Space Sci. Rev., 214, 52 [NASA ADS] [CrossRef] [Google Scholar]
  8. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
  9. Blum, J., Gundlach, B., Mühle, S., & Trigo-Rodriguez, J. M. 2014, Icarus, 235, 156 [Google Scholar]
  10. Blum, J., Gundlach, B., Krause, M., et al. 2017, MNRAS, 469, S755 [Google Scholar]
  11. Blum, J., Bischoff, D., & Gundlach, B. 2022, Universe, 8, 381 [CrossRef] [Google Scholar]
  12. Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 175, 111 [Google Scholar]
  13. Britt, D. T., Boice, D. C., Buratti, B. J., et al. 2004, Icarus, 167, 45 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brown, M. E. 2001, AJ, 121, 2804 [Google Scholar]
  15. Buratti, B. J., Hicks, M. D., Soderblom, L. A., et al. 2004, Icarus, 167, 16 [NASA ADS] [CrossRef] [Google Scholar]
  16. Desch, S. J., & Jackson, A. P. 2021, J. Geophys. Res. Planets, 126, e06807 [NASA ADS] [CrossRef] [Google Scholar]
  17. Desch, S. J., & Jackson, A. P. 2022, Astrobiology, 22, 1400 [NASA ADS] [CrossRef] [Google Scholar]
  18. Duxbury, T. C., Newburn, R. L., & Brownlee, D. E. 2004, J. Geophys. Res. Planets, 109, E12S02 [Google Scholar]
  19. Farnham, T. L., & Thomas, P. C. 2013, Plate Shape Model of Comet 103P/HARTLEY 2 V1.0 [Google Scholar]
  20. Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2023, Planet. Sci. J., 4, 29 [NASA ADS] [CrossRef] [Google Scholar]
  21. Flekkøy, E. G., Luu, J., & Toussaint, R. 2019, ApJ, 885, L41 [CrossRef] [Google Scholar]
  22. Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, Nat. Astron., 1, 0088 [Google Scholar]
  23. Frezzotti, A. 1998, Comput. Math. Applic., 35, 103 [CrossRef] [Google Scholar]
  24. Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216 [Google Scholar]
  25. Gomes, R., Nesvorný, D., Morbidelli, A., Deienno, R., & Nogueira, E. 2018, Icarus, 306, 319 [NASA ADS] [CrossRef] [Google Scholar]
  26. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [Google Scholar]
  27. Harmon, J. K., Nolan, M. C., Giorgini, J. D., & Howell, E. S. 2010, Icarus, 207, 499 [CrossRef] [Google Scholar]
  28. Hong, A., & Morris, A. 2021, Granular Matter, 23, 71 [CrossRef] [Google Scholar]
  29. Jackson, A. P., & Desch, S. J. 2021, J. Geophys. Res. Planets, 126, e06706 [NASA ADS] [CrossRef] [Google Scholar]
  30. Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36 [NASA ADS] [CrossRef] [Google Scholar]
  31. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  32. Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, 547 [Google Scholar]
  33. Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257 [Google Scholar]
  34. Jutzi, M., & Asphaug, E. 2015, Science, 348, 1355 [NASA ADS] [CrossRef] [Google Scholar]
  35. Jutzi, M., & Benz, W. 2017, A&A, 597, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Kavelaars, J. J., Petit, J.-M., Gladman, B., et al. 2021, ApJ, 920, L28 [NASA ADS] [CrossRef] [Google Scholar]
  37. Keane, J. T., Porter, S. B., Beyer, R. A., et al. 2022, J. Geophys. Res. Planets, 127, e07068 [NASA ADS] [CrossRef] [Google Scholar]
  38. Keller, H. U., Arpigny, C., Barbieri, C., et al. 1986, Nature, 321, 320 [NASA ADS] [CrossRef] [Google Scholar]
  39. Keller, H. U., Delamere, W. A., Reitsema, H. J., Huebner, W. F., & Schmidt, H. U. 1987, A&A, 187, 807 [NASA ADS] [Google Scholar]
  40. Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Lacerda, P. 2011, AJ, 142, 90 [NASA ADS] [CrossRef] [Google Scholar]
  42. Lacerda, P., McNeill, A., & Peixinho, N. 2014, MNRAS, 437, 3824 [NASA ADS] [CrossRef] [Google Scholar]
  43. Lorek, S., Lacerda, P., & Blum, J. 2018, A&A, 611, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Ma, D., & Ahmadi, G. 1986, J. Chem. Phys., 84, 3449 [CrossRef] [Google Scholar]
  45. Marohnic, J. C., Richardson, D. C., McKinnon, W. B., et al. 2021, Icarus, 356, 113824 [NASA ADS] [CrossRef] [Google Scholar]
  46. Mashchenko, S. 2019, MNRAS, 489, 3003 [NASA ADS] [CrossRef] [Google Scholar]
  47. Massironi, M., Simioni, E., Marzari, F., et al. 2015, Nature, 526, 402 [NASA ADS] [CrossRef] [Google Scholar]
  48. McKinnon, W. B., Richardson, D. C., Marohnic, J. C., et al. 2020, Science, 367, aay6620 [NASA ADS] [CrossRef] [Google Scholar]
  49. Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378 [Google Scholar]
  50. Merenyl, E., Foldy, L., Szego, K., Toth, I., & Kondor, A. 1990, Icarus, 86, 9 [NASA ADS] [CrossRef] [Google Scholar]
  51. Michel, P., & Richardson, D. C. 2013, A&A, 554, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Michel, P., Ballouz, R. L., Barnouin, O. S., et al. 2020, Nat. Commun., 11, 2655 [NASA ADS] [CrossRef] [Google Scholar]
  53. Micheli, M., Farnocchia, D., Meech, K. J., et al. 2018, Nature, 559, 223 [Google Scholar]
  54. Montanero, J. M., & Santos, A. 1996, Phys. Rev. E, 54, 438 [NASA ADS] [CrossRef] [Google Scholar]
  55. Montanero, J. M., & Santos, A. 1997, Phys. Fluids, 9, 2057 [NASA ADS] [CrossRef] [Google Scholar]
  56. Morbidelli, A., & Rickman, H. 2015, A&A, 583, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785 [Google Scholar]
  58. Nesvorný, D., Vokrouhlický, D., Bottke, W. F., Noll, K., & Levison, H. F. 2011, AJ, 141, 159 [Google Scholar]
  59. Nesvorný, D., Parker, J., & Vokrouhlický, D. 2018, AJ, 155, 246 [CrossRef] [Google Scholar]
  60. Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nat. Astron., 3, 808 [CrossRef] [Google Scholar]
  61. Nesvorný, D., Li, R., Simon, J. B., et al. 2021, Planet. Sci. J., 2, 27 [CrossRef] [Google Scholar]
  62. Noll, K. S., Grundy, W. M., Chiang, E. I., Margot, J. L., & Kern, S. D. 2008a, in The Solar System Beyond Neptune, eds. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson, 345 [Google Scholar]
  63. Noll, K. S., Grundy, W. M., Stephens, D. C., Levison, H. F., & Kern, S. D. 2008b, Icarus, 194, 758 [NASA ADS] [CrossRef] [Google Scholar]
  64. Oberst, J., Giese, B., Howington-Kraus, E., et al. 2004, Icarus, 167, 70 [NASA ADS] [CrossRef] [Google Scholar]
  65. Onoda, G. Y., & Liniger, E. G. 1990, Phys. Rev. Lett., 64, 2727 [CrossRef] [PubMed] [Google Scholar]
  66. ‘Oumuamua ISSI Team, Bannister, M. T., Bhandare, A., et al. 2019, Nat. Astron., 3, 594 [NASA ADS] [CrossRef] [Google Scholar]
  67. Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131 [NASA ADS] [CrossRef] [Google Scholar]
  68. Polak, B., & Klahr, H. 2023, ApJ, 943, 125 [NASA ADS] [CrossRef] [Google Scholar]
  69. Pöschel, T., & Schwager, T. 2005, Computational Granular Dynamics (Springer) [Google Scholar]
  70. Rafikov, R. R. 2018, ApJ, 867, L17 [NASA ADS] [CrossRef] [Google Scholar]
  71. Raymond, S. N., Armitage, P. J., & Veras, D. 2018, ApJ, 856, L7 [Google Scholar]
  72. Robinson, J. E., Fraser, W. C., Fitzsimmons, A., & Lacerda, P. 2020, A&A, 643, A55 [EDP Sciences] [Google Scholar]
  73. Safrit, T. K., Steckloff, J. K., Bosh, A. S., et al. 2021, Planet. Sci. J., 2, 14 [NASA ADS] [CrossRef] [Google Scholar]
  74. Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69 [Google Scholar]
  75. Schwartz, S. R., Michel, P., Jutzi, M., et al. 2018, Nat. Astron., 2, 379 [NASA ADS] [CrossRef] [Google Scholar]
  76. Seligman, D., & Laughlin, G. 2020, ApJ, 896, L8 [NASA ADS] [CrossRef] [Google Scholar]
  77. Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, Planet. Sci. J., 4, 35 [NASA ADS] [CrossRef] [Google Scholar]
  78. Sheppard, S. S., & Jewitt, D. 2004, AJ, 127, 3023 [NASA ADS] [CrossRef] [Google Scholar]
  79. Showalter, M. R., Benecchi, S. D., Buie, M. W., et al. 2021, Icarus, 356, 114098 [NASA ADS] [CrossRef] [Google Scholar]
  80. Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044 [Google Scholar]
  81. Simon, J. B., Blum, J., Birnstiel, T., & Nesvorný, D. 2022, ArXiv e-prints [arXiv:2212.04509] [Google Scholar]
  82. Stern, S. A., Weaver, H. A., Spencer, J. R., et al. 2019, Science, 364, aaw9771 [NASA ADS] [CrossRef] [Google Scholar]
  83. Thirouin, A., & Sheppard, S. S. 2018, AJ, 155, 248 [NASA ADS] [CrossRef] [Google Scholar]
  84. Thirouin, A., & Sheppard, S. S. 2019, AJ, 157, 228 [NASA ADS] [CrossRef] [Google Scholar]
  85. Thirouin, A., Sheppard, S. S., & Noll, K. S. 2017, ApJ, 844, 135 [NASA ADS] [CrossRef] [Google Scholar]
  86. Thomas, P. C., A’Hearn, M. F., Veverka, J., et al. 2013, Icarus, 222, 550 [NASA ADS] [CrossRef] [Google Scholar]
  87. Torquato, S., & Stillinger, F. H. 2007, J. Appl. Phys., 102 [CrossRef] [Google Scholar]
  88. Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261 [Google Scholar]
  89. Van Beijeren, H., & Ernst, M. H. 1973, Physica, 68, 437 [NASA ADS] [CrossRef] [Google Scholar]
  90. Vavilov, D. E., Eggl, S., Medvedev, Y. D., & Zatitskiy, P. B. 2019, A&A, 622, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Wahlberg Jansson, K., & Johansen, A. 2014, A&A, 570, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Weidling, R., & Blum, J. 2015, Icarus, 253, 31 [NASA ADS] [CrossRef] [Google Scholar]
  93. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
  94. Zhao, Y., Rezac, L., Skorov, Y., et al. 2021, Nat. Astron., 5, 139 [Google Scholar]
  95. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1

Comet and KBO sizes.

In the text
Table D.1

Ellipsoidal fit 103P/Hartley 2.

In the text

All Figures

thumbnail Fig. 1

Shapes of planetesimals for different values of the initial angular momentum. The rows show the column densities of cells that have a volume-filling factor ≥0.5 in the xy- and xz-planes for different values of the initial angular momentum of the cloud at the end of each simulation. We omit the yz-plane because of the evident cylindrical symmetry. The contour of the ellipsoid fitted to the planetesimals in the respective plane are shown in black with the centre of the ellipse marked with a black circle. The figure shows the results for initial values of L/LJ = 1.5, 1.0, 0.5 and 0.0 to illustrate the effect of angular momentum of the cloud on the flattening of the planetesimals.
In the text
thumbnail Fig. 2

Axes ratios of planetesimals as a function of the initial angular momentum of the cloud. The axes ratios of the two lobes of Arrokoth are shown as dashed (big lobe) and dotted (small lobe) lines.
In the text
thumbnail Fig. 3

Average axes ratios of planetesimals for the initial angular momentum of the cloud in the range 0.5 to 0.9. For each initial L/LJ. the results of 5 different realisations are averaged. The axes ratios of the two lobes of Arrokoth are shown as dashed (big lobe) and dotted (small lobe) lines.
In the text
thumbnail Fig. 4

Angular momentum of planetesimals as a function of the initial angular momentum of the cloud. The dashed line indicates where the final angular momentum of the planetesimal Lfinal equals the initial angular momentum and the dotted line where Lfinal equals half of the initial value. The colour coding corresponds to the mass that ends up in the planetesimal. The red crosses show the total angular momentum Llot in the simulation.
In the text
thumbnail Fig. 5

Axes ratios of Solar System bodies, ‘Oumuamua, and simulated planetesimals. The simulated planetesimals are shown as circles with a colour that indicates the initial angular momentum of the pebble cloud L/LJ. For L/LJ > 1, we used open symbols. The data points with error bars are from simulations where we averaged over five runs with varying random seed. The Solar System bodies and ‘Oumuamua are shown with open diamond symbols and different colours: 1P/Halley (blue), 8P/Tuttle (orange), 9P/Tempel 1 (green), 19P/Borrelly (red), 67P/Churyumov-Gerasimenko (purple), 81P/Wild 2 (brown), 103P/Hartley 2 (pink), Arrokoth (grey), 1I/2017 U1 (‘Oumuamua) (olive). Left: axes ratios of the whole bodies. Right: axes ratios of the individual lobes of the bilobed and contact-binary comets, Arrokoth, and the simulated planetesimals. We plot the axes ratios of the individual lobes with an upper triangle for the big lobe and a lower triangle for the small lobe. For Hartley 2, we included the two lobes that we found from analysing the shape model (see Appendix D).
In the text
thumbnail Fig. E.1

Shape model of a planetesimal with L/LJ=0.5. Each volume element with a volume-filling factor ≥0.5 is plotted to visualise the shape of the planetesimal. The centres of the volume elements are projected into the xy-, xz-, and yz-planes to visualise the shape of the planetesimal that can not be seen in the three-dimensional representation.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.