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/00046361/202347742  
Published online  04 March 2024 
Formation of flattened planetesimals by gravitational collapse of rotating pebble clouds
^{1}
Centre for Star and Planet Formation, Globe Institute, University of Copenhagen,
Øster Voldgade 5–7,
1350
Copenhagen,
Denmark
email: sebastian.lorek@sund.ku.dk
^{2}
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43,
221 00
Lund,
Sweden
Received:
17
August
2023
Accepted:
13
December
2023
Planetesimals are believed to form by the gravitational collapse of aerodynamically concentrated clumps of pebbles. Many properties of the objects in the cold classical Kuiper belt – such as binarity, rotation, and size distribution – are in agreement with this gravitational collapse model. Further support comes from the pebblepile structure inferred for comet nuclei. For this study, we simulated the final assembly of a planetesimal from the gravitational collapse of a rotating clump of pebbles. We implemented a numerical method from granular dynamics to follow the collapse that includes the transition from a pebble swarm to solid cells at a high density. We compared the shapes of the simulated planetesimals with the shapes of the lobes of contact binaries and bilobed Solar System objects. We find that the gravitational collapse of slowly rotating pebble clouds naturally explains the formation of flattened ellipsoidal bodies. This result agrees well with the flattened structure of the bilobed planetesimal Arrokoth and the shapes of the components of bilobed comets.
Key words: methods: numerical / planets and satellites: formation
© The Authors 2024
Open 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 decimetresized aggregates that have grown by collisions of micrometresized 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 centimetresize 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 filamentlike structures that subsequently fragment under their own gravity to form pebblepile 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 selfgravity 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 byproducts 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 reaccumulation of the fragments from catastrophic collisions of larger bodies (Michel & Richardson 2013; Michel et al. 2020). Therefore, the irregular shapes of kilometresized 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 activitydriven erosion changes the surface morphology of the nuclei. Anisotropic mass loss caused by nonuniform insolation of the nucleus is capable of producing irregularshaped 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 cometsized planetesimals during the dynamical evolution of the outer Solar System suggests that kilometresized 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 reaccumulation of the fragments produced in subcatastrophic 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 socalled cold classical Kuiper belt, forms a narrow ring between approximately 42.4 au and 47.7 au, with a subpopulation centred at ~44 au (the socalled 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 initialmass 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 highresolution 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 QG_{298}, a contactbinary fraction of ≳ 10–30% was estimated (Sheppard & Jewitt 2004; Lacerda 2011). Other candidates are the Kuiper belt objects 2003 SQ_{317} and 2004 TT_{357} 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 contactbinary 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 reaccumulation, comets and the cold classical Kuiper belt objects might still be representative of the shapes of the kilometresized primordial planetesimals that formed through the gravitational collapse of pebble clouds. The dynamics of the pebble cloud collapse has been studied in high resolution Nbody simulations or by means of hydrodynamic simulations of the selfgravitating 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 kilometresized planetesimal from the gravitational collapse of a rotating clump of pebbles using a MonteCarlo method for granular dynamics and show that, depending on the angular momentum content, the final shapes of kilometresized planetesimals are in agreement with the shapes of the leftover 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 MonteCarlo 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 wellstudied 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 × 10^{14} 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 massfraction of 75% with an average rock density of 3500 kg m^{−3}), the volumefilling 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 lowvelocity 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 irregularshaped kilometresized 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 Jupiterfamily comet 67P/ChuryumovGerasimenko for almost 2 yr along its orbit, obtained highresolution images and detailed shape information of a kilometresized 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 × 10^{12} 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 bestfit model of Mashchenko (2019) suggests a thin disc with dimensions 115 × 111 × 19 m. A cigarshaped 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 lightcurve 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 kilometresized planetesimal requires ~10^{18} millimetresized 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 selfgravity, or in other words, solve the Boltzmann equation. We do this here by using a directsimulation 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). Selfgravity is accounted for in the advection step by using a driftkickdrift 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 twoparticle distribution function can be written as the product of two oneparticle distribution functions (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 twoparticle correlation function, such that the distribution function reads (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 hardsphere fluid (Ma & Ahmadi 1986; see Appendix A). The Boltzmann equation then turns into the BoltzmannEnskog 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 overpacking of cells and further allows to identify fully packed cells as the solid elements comprising the planetesimal. The maximum volumefilling 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 meanfree 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 kilometresized planetesimal from a pebble cloud with low angular momentum content. We cannot resolve the entire collapse from a Hillradiussized pebble cloud down to a planetesimal because the required number of grid cells to achieve subkilometre resolution would be too high. Instead, we assume that the pebble cloud has already fragmented into subclumps that would eventually form, for example, the components of a binary as has been shown in Nbody simulations of the gravitational collapse (e.g. Polak & Klahr 2023).
Fig. 1 Shapes of planetesimals for different values of the initial angular momentum. The rows show the column densities of cells that have a volumefilling factor ≥0.5 in the xy and xzplanes for different values of the initial angular momentum of the cloud at the end of each simulation. We omit the yzplane 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/L_{J} = 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 R_{p} = 1 km (as chosen in this study) this meant that the initial cloud radius was R_{c} = 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 , 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 M_{p} and effective radius R_{p} has an angular momentum of (3)
We therefore set the angular momentum content in units of L_{J} because the cloud can collapse into a single body only for L/L_{J} ≲ 1. To study the effect of angular momentum on the planetesimal shape we systematically varied L/L_{J} between 0 (no rotation) and 2 (supercritical rotation).
We conducted the study for a pebble cloud that had the mass of a planetesimal of R_{p} = 1 km and density 1280kg m^{−3}, which corresponded to a volumefilling factor of the planetesimal of 0.64 and a material density of the pebbles of 2000 kg m^{3}. 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 freefall 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 freefall times of the cloud until a solid body formed that no longer changed its shape significantly. The freefall time was of the order of a day. We then identified all cells that had a volumefilling 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 volumefilling 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 zerogravity limit with volumefilling 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 zdirection which is parallel to the initial rotation axis. In the plane perpendicular to the rotation axis (xyplane), the planetesimals remain more round. As expected for very high angular momentum content L/L_{J} ≳ 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/L_{J} ≲ 1.5, very thin disclike objects form. When the angular momentum is L/L_{J} ≲ 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 disclike and lenticularshaped planetesimals in our simulations. However, we do not see elongated ellipsoids that would resemble bodies such as 19P/Borrelly or 103P/Hartley 2.
Comet and KBO sizes.
4.3 Axes ratios of planetesimals
We used linear leastsquares to fit an ellipsoid that is centred at the centreofmass 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 baplane and flattened in the ca and cbplanes 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/L_{J} ≲ 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/L_{J} ≲ 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/L_{J} in the range 0.5–0.9. For each value of L/L_{J}, we averaged the axes ratios over five realisations for the same initial conditions to obtain a measure for the uncertainty of our MonteCarlo simulations. Figure 3 shows that the uncertainty of the axes ratios are ~10%.
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. 
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/L_{J}. 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. 
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 L_{final} equals the initial angular momentum and the dotted line where L_{final} 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 L_{lot} 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 kilometresized 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 disclike shape of ‘Oumuamua would be consistent with the very flat planetesimals that form in pebble clouds with 1 ≲ L/L_{J} ≲ 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/CG, 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 Nbody simulations and shows that the pebble cloud typically fragments into a number of subclumps 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 lowvelocity 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 subclumps. 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.
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/L_{J}. For L/L_{J} > 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/ChuryumovGerasimenko (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 contactbinary 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 nongravitational acceleration (Micheli et al. 2018). It is unlikely that the nongravitational acceleration is caused by strong volatile outgassing, 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 H_{2} could explain the nongravitational acceleration. Desch & Jackson (2021, 2022) and Jackson & Desch (2021) argue that ‘Oumuamua could be a nitrogenice fragment excavated from an exoPluto with its shape heavily affected by N_{2}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 cometlike planetesimal. They show that the nongravitational 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 subkilometresized nearEarth asteroids that exhibit significant nongravitational 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 disclike 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 kilometresized planetesimal, we cannot follow the entire collapse from a Hillsized 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 subclumps that form during the collapse (Nesvorný et al. 2021; Polak & Klahr 2023). Because of the fixed grid, cells that are filled to the maximum volumefilling 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 disclike flat structures (L/L_{J} ≳ 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 angularmomentum clouds that collapse into a single body.
6 Conclusion
We developed a MonteCarlo 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/L_{J} ≲ 1, where L_{J} 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/L_{J} ≲ 1. Lastly, the disclike 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/L_{J} ≲ 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 724687PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 201804867), 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 N_{swm} = N_{phys}/N pebbles, where N_{phys} =M_{p}/m_{peb} is the total number of pebbles of mass m_{peb} that is needed for a planetesimal of mass M_{p}. Each particle has position x, velocity v, mass m, and radius R. The MonteCarlo 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 (A.1)
Here, σ is the diameter of a pebble in case of a monodisperse size distribution of pebbles; for a nonuniform size distribution, σ is chosen at random in the interval between the minimum and the maximum pebble radius such that R_{i} + R_{min}≤σ≤R_{i}+R_{max}. The relative velocity between particles i and j is v_{ij} = v_{j} − v_{i} and e_{i} is the direction of the collision, that is, the unit vector pointing from particle i to j. Finally, χ(x_{i}, x_{i} + σe_{i}) is the local equilibrium pair correlation function (or Enskog factor) (Van Beijeren & Ernst 1973), which depends on the density in the cell J, and n_{J} 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 hardsphere equation of state from Ma & Ahmadi (1986), which gives (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 (A.3)
where ω_{max}= max_{ij}({ω_{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 e_{i} is chosen at random with equiprobability.
 3.
A particle j from the cell J where x_{i}+σe_{i} points to is chosen at random.
 4.
The collision between particles i and j is accepted with probability Θ(e_{i} ⋅ v_{ij})ω_{ij}/ω_{max}. Here, Θ(e_{i} ⋅ v_{ij}) 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 , 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 Selfgravity of the cloud
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 (B.2) (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 (B.4)
We first calculate using the Fast Fourier Transform (FFT), then solve the above equation, and in a last step calculate the inverse FFT of 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 finitedifference scheme (B.5) (B.6) (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 (C.1)
We have N data points (x_{i}, y_{i}, z_{i}) 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 , we can write the matrix equation (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 righthand 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 (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, (C.4)
where we introduced the vector y=(x y z 1)^{⊤} and the matrix Q (C.5)
To find the centre x_{0} = (x_{0} y_{0} z_{0}) of the ellipsoid, we introduce a translation matrix of the form (C.6)
such that Ty moves the ellipsoid to the centre x_{0}. 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 (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 T_{14} to T_{34} arise from the fact that the ellipsoid is not centred at (0, 0, 0) but shifted to x_{0}, 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 T_{44} is the normalisation for recovering the enforced J = −1. The individual entries are (C.8) (C.9) (C.10) (C.11)
To find the centre of the ellipsoid, we notice that the entries T_{14}=T_{24}=T_{34}=0 need to be equal to zero. This gives the following linear matrix equation (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 (C.13)
which we need to normalise by dividing all entries by −T_{44} in order to recover J = −1. The eigenvalues are the solution of the characteristic polynomial (C.14)
where 𝕀_{3} is the 3 × 3 identity matrix and λ is the eigenvalue. Because in diagonal form the matrix describing the ellipsoid is (C.15)
which results in λ_{a}x^{2} + λ_{b}y^{2} + λ_{c}z^{2} = 1, the axes are related to the eigenvalues as , and . 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.
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/L_{J} = 0.5. One can see the flattened structure of the planetesimal arising from angular momentum conservation.
Fig. E.1 Shape model of a planetesimal with L/L_{J}=0.5. Each volume element with a volumefilling 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 yzplanes to visualise the shape of the planetesimal that can not be seen in the threedimensional representation. 
References
 A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258 [Google Scholar]
 A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396 [CrossRef] [Google Scholar]
 Alexander, F. J., & Garcia, A. L. 1997, Comput. Phys., 11, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Brown, M. E., & Fraser, W. C. 2011, ApJ, 738, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Bergner, J. B., & Seligman, D. Z. 2023, Nature, 615, 610 [NASA ADS] [CrossRef] [Google Scholar]
 Bird, G. A. 1994, Molecular Gas Dynamics And The Direct Simulation Of Gas Flows (Oxford: Clarendon Press) [Google Scholar]
 Blum, J. 2018, Space Sci. Rev., 214, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2014, Icarus, 235, 156 [Google Scholar]
 Blum, J., Gundlach, B., Krause, M., et al. 2017, MNRAS, 469, S755 [Google Scholar]
 Blum, J., Bischoff, D., & Gundlach, B. 2022, Universe, 8, 381 [CrossRef] [Google Scholar]
 Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 175, 111 [Google Scholar]
 Britt, D. T., Boice, D. C., Buratti, B. J., et al. 2004, Icarus, 167, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E. 2001, AJ, 121, 2804 [Google Scholar]
 Buratti, B. J., Hicks, M. D., Soderblom, L. A., et al. 2004, Icarus, 167, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Desch, S. J., & Jackson, A. P. 2021, J. Geophys. Res. Planets, 126, e06807 [NASA ADS] [CrossRef] [Google Scholar]
 Desch, S. J., & Jackson, A. P. 2022, Astrobiology, 22, 1400 [NASA ADS] [CrossRef] [Google Scholar]
 Duxbury, T. C., Newburn, R. L., & Brownlee, D. E. 2004, J. Geophys. Res. Planets, 109, E12S02 [Google Scholar]
 Farnham, T. L., & Thomas, P. C. 2013, Plate Shape Model of Comet 103P/HARTLEY 2 V1.0 [Google Scholar]
 Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2023, Planet. Sci. J., 4, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Flekkøy, E. G., Luu, J., & Toussaint, R. 2019, ApJ, 885, L41 [CrossRef] [Google Scholar]
 Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, Nat. Astron., 1, 0088 [Google Scholar]
 Frezzotti, A. 1998, Comput. Math. Applic., 35, 103 [CrossRef] [Google Scholar]
 Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216 [Google Scholar]
 Gomes, R., Nesvorný, D., Morbidelli, A., Deienno, R., & Nogueira, E. 2018, Icarus, 306, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [Google Scholar]
 Harmon, J. K., Nolan, M. C., Giorgini, J. D., & Howell, E. S. 2010, Icarus, 207, 499 [CrossRef] [Google Scholar]
 Hong, A., & Morris, A. 2021, Granular Matter, 23, 71 [CrossRef] [Google Scholar]
 Jackson, A. P., & Desch, S. J. 2021, J. Geophys. Res. Planets, 126, e06706 [NASA ADS] [CrossRef] [Google Scholar]
 Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, 547 [Google Scholar]
 Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257 [Google Scholar]
 Jutzi, M., & Asphaug, E. 2015, Science, 348, 1355 [NASA ADS] [CrossRef] [Google Scholar]
 Jutzi, M., & Benz, W. 2017, A&A, 597, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kavelaars, J. J., Petit, J.M., Gladman, B., et al. 2021, ApJ, 920, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Keane, J. T., Porter, S. B., Beyer, R. A., et al. 2022, J. Geophys. Res. Planets, 127, e07068 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Arpigny, C., Barbieri, C., et al. 1986, Nature, 321, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, H. U., Delamere, W. A., Reitsema, H. J., Huebner, W. F., & Schmidt, H. U. 1987, A&A, 187, 807 [NASA ADS] [Google Scholar]
 Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lacerda, P. 2011, AJ, 142, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Lacerda, P., McNeill, A., & Peixinho, N. 2014, MNRAS, 437, 3824 [NASA ADS] [CrossRef] [Google Scholar]
 Lorek, S., Lacerda, P., & Blum, J. 2018, A&A, 611, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ma, D., & Ahmadi, G. 1986, J. Chem. Phys., 84, 3449 [CrossRef] [Google Scholar]
 Marohnic, J. C., Richardson, D. C., McKinnon, W. B., et al. 2021, Icarus, 356, 113824 [NASA ADS] [CrossRef] [Google Scholar]
 Mashchenko, S. 2019, MNRAS, 489, 3003 [NASA ADS] [CrossRef] [Google Scholar]
 Massironi, M., Simioni, E., Marzari, F., et al. 2015, Nature, 526, 402 [NASA ADS] [CrossRef] [Google Scholar]
 McKinnon, W. B., Richardson, D. C., Marohnic, J. C., et al. 2020, Science, 367, aay6620 [NASA ADS] [CrossRef] [Google Scholar]
 Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378 [Google Scholar]
 Merenyl, E., Foldy, L., Szego, K., Toth, I., & Kondor, A. 1990, Icarus, 86, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, P., & Richardson, D. C. 2013, A&A, 554, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michel, P., Ballouz, R. L., Barnouin, O. S., et al. 2020, Nat. Commun., 11, 2655 [NASA ADS] [CrossRef] [Google Scholar]
 Micheli, M., Farnocchia, D., Meech, K. J., et al. 2018, Nature, 559, 223 [Google Scholar]
 Montanero, J. M., & Santos, A. 1996, Phys. Rev. E, 54, 438 [NASA ADS] [CrossRef] [Google Scholar]
 Montanero, J. M., & Santos, A. 1997, Phys. Fluids, 9, 2057 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., & Rickman, H. 2015, A&A, 583, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785 [Google Scholar]
 Nesvorný, D., Vokrouhlický, D., Bottke, W. F., Noll, K., & Levison, H. F. 2011, AJ, 141, 159 [Google Scholar]
 Nesvorný, D., Parker, J., & Vokrouhlický, D. 2018, AJ, 155, 246 [CrossRef] [Google Scholar]
 Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nat. Astron., 3, 808 [CrossRef] [Google Scholar]
 Nesvorný, D., Li, R., Simon, J. B., et al. 2021, Planet. Sci. J., 2, 27 [CrossRef] [Google Scholar]
 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]
 Noll, K. S., Grundy, W. M., Stephens, D. C., Levison, H. F., & Kern, S. D. 2008b, Icarus, 194, 758 [NASA ADS] [CrossRef] [Google Scholar]
 Oberst, J., Giese, B., HowingtonKraus, E., et al. 2004, Icarus, 167, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Onoda, G. Y., & Liniger, E. G. 1990, Phys. Rev. Lett., 64, 2727 [CrossRef] [PubMed] [Google Scholar]
 ‘Oumuamua ISSI Team, Bannister, M. T., Bhandare, A., et al. 2019, Nat. Astron., 3, 594 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Polak, B., & Klahr, H. 2023, ApJ, 943, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Pöschel, T., & Schwager, T. 2005, Computational Granular Dynamics (Springer) [Google Scholar]
 Rafikov, R. R. 2018, ApJ, 867, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Armitage, P. J., & Veras, D. 2018, ApJ, 856, L7 [Google Scholar]
 Robinson, J. E., Fraser, W. C., Fitzsimmons, A., & Lacerda, P. 2020, A&A, 643, A55 [EDP Sciences] [Google Scholar]
 Safrit, T. K., Steckloff, J. K., Bosh, A. S., et al. 2021, Planet. Sci. J., 2, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Schäfer, U., Yang, C.C., & Johansen, A. 2017, A&A, 597, A69 [Google Scholar]
 Schwartz, S. R., Michel, P., Jutzi, M., et al. 2018, Nat. Astron., 2, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Seligman, D., & Laughlin, G. 2020, ApJ, 896, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, Planet. Sci. J., 4, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Sheppard, S. S., & Jewitt, D. 2004, AJ, 127, 3023 [NASA ADS] [CrossRef] [Google Scholar]
 Showalter, M. R., Benecchi, S. D., Buie, M. W., et al. 2021, Icarus, 356, 114098 [NASA ADS] [CrossRef] [Google Scholar]
 Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044 [Google Scholar]
 Simon, J. B., Blum, J., Birnstiel, T., & Nesvorný, D. 2022, ArXiv eprints [arXiv:2212.04509] [Google Scholar]
 Stern, S. A., Weaver, H. A., Spencer, J. R., et al. 2019, Science, 364, aaw9771 [NASA ADS] [CrossRef] [Google Scholar]
 Thirouin, A., & Sheppard, S. S. 2018, AJ, 155, 248 [NASA ADS] [CrossRef] [Google Scholar]
 Thirouin, A., & Sheppard, S. S. 2019, AJ, 157, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Thirouin, A., Sheppard, S. S., & Noll, K. S. 2017, ApJ, 844, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, P. C., A’Hearn, M. F., Veverka, J., et al. 2013, Icarus, 222, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Torquato, S., & Stillinger, F. H. 2007, J. Appl. Phys., 102 [CrossRef] [Google Scholar]
 Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261 [Google Scholar]
 Van Beijeren, H., & Ernst, M. H. 1973, Physica, 68, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Vavilov, D. E., Eggl, S., Medvedev, Y. D., & Zatitskiy, P. B. 2019, A&A, 622, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wahlberg Jansson, K., & Johansen, A. 2014, A&A, 570, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidling, R., & Blum, J. 2015, Icarus, 253, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
 Zhao, Y., Rezac, L., Skorov, Y., et al. 2021, Nat. Astron., 5, 139 [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Shapes of planetesimals for different values of the initial angular momentum. The rows show the column densities of cells that have a volumefilling factor ≥0.5 in the xy and xzplanes for different values of the initial angular momentum of the cloud at the end of each simulation. We omit the yzplane 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/L_{J} = 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 
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 
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/L_{J}. 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 
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 L_{final} equals the initial angular momentum and the dotted line where L_{final} 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 L_{lot} in the simulation. 

In the text 
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/L_{J}. For L/L_{J} > 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/ChuryumovGerasimenko (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 contactbinary 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 
Fig. E.1 Shape model of a planetesimal with L/L_{J}=0.5. Each volume element with a volumefilling 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 yzplanes to visualise the shape of the planetesimal that can not be seen in the threedimensional representation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.