Formation of flattened planetesimals by gravitational collapse of rotating pebble clouds

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 pebble-pile 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.


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. 2007Johansen et al. , 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 primor-dial, 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(Blum et al. , 2017(Blum et al. , 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. 2005A'Hearn et al. , 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 reaccumulation 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(Nesvorný et al. , 2019;;Robinson et al. 2020;Nesvorný et al. 2021;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 contact-binary 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 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(Nesvorný et al. , 2019;;Robinson et al. 2020;Nesvorný et al. 2021;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 kilometresized planetesimals are in agreement with the shapes of the leftover planetesimals in the Solar System.
The paper is organised as follows.In Section 2, we broadly summarise available shape information of minor bodies.In Section 3, we provide an overview of the Monte-Carlo method, with a more detailed description in the Appendix A. In Section 4, we go throught the results of our simulations.In Section 5, we discuss our results in the context of the Solar System bodies.Finally, in Section 6, we summarise the main findings of our study.

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 1I/2017 U1 ('Oumuamua).

Arrokoth
The overall dimensions of Arrokoth are 35.95×19.90×9.75km.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.05km for the large lobe Wenu and 15.75×13.85×9.75km 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 Notes.Axis dimensions a, b, and c are diameters.Fit refers to how the axis dimensions were estimated: shape model (shape), triaxial ellipsoid fit (ellipsoid), or prolate ellipsoid fit (prolate).For prolate models, the two small axes b and c are equal.A shape model is a 3D reconstruction of the shape of the object based on images taken by the spacecraft.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.

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(Keller et al. , 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. 2005A'Hearn et al. , 2011)).Lastly, the most recent mission to a comet, ESA's Rosetta mission, which accompanied the Jupiter-family comet 67P/Churyumov-Gerasimenko for almost two years 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×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 obser-vations of comet 8P/Tuttle revealed a contact binary with a large lobe of size 5.75×4.11km and a small lobe of size 4.25×3.27km (Harmon et al. 2010).

1I/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.

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 ∼10 18 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 (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 (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.

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).

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 √ k B T/m≈(1/3)v vir , where v 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 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 (super-critical rotation).
We conducted the study for a pebble cloud that had the mass of a planetesimal of R p =1 km and density 1280 kg 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 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 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 freefall 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).

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/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. 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 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.

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/L J ≲1 the planetesimals become more spherical as is already indicated in Figure 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 Monte-Carlo simulations.Figure 3 shows that the uncertainty of the axes ratios are ∼10 %.

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.

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 Section 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.

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/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/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 1P/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(Blum et al. , 2017(Blum et al. , 2022;;Nesvorný et al. 2010Nesvorný et al. , 2019Nesvorný et al. , 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;Robinson et al. 2020;Nesvorný et al. 2021;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.

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 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 non-gravitational acceleration.Desch & Jackson (2021), Jackson &Desch (2021), andDesch &Jackson (2022) argue that 'Oumuamua could be a nitrogen-ice fragment excavated from an exo-Pluto 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 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 subkilometre-sized near-Earth 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 disc-like flattened shape of 'Oumuamua could arise naturally from the gravitational collapse of a pebble cloud with high angular momentum.

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 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/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 angular-momentum clouds that collapse into a single body.

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/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 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/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.
Sebastian Lorek, Anders Johansen: Formation of flattened planetesimals by gravitational collapse of rotating pebble clouds 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 (x 2 i y 2 i z 2 i x i y i x i z i y i z i x i y i z i ), we can write the matrix equation 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, (C.4) where we introduced the vector y=(x y z 1) ⊤ and the matrix Q To find the centre x 0 =(x 0 y 0 z 0 ) of the ellipsoid, we introduce a translation matrix of the form 1 0 0 x 0 0 1 0 y 0 0 0 1 z 0 0 0 0 1 such that T y 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 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.T 44 = Ax 2 0 + By 2 0 + Cz 2 0 + Dx 0 y 0 + Ex 0 z 0 + Fy 0 z 0 + Gx 0 + Hy 0 + Iz 0 + J. (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 where I 3 is the 3×3 identity matrix and λ is the eigenvalue.Because in diagonal form the matrix describing the ellipsoid is which results in λ a x 2 +λ b y 2 +λ c z 2 =1, the axes are related to the eigenvalues as a=λ −1/2 a , b=λ −1/2 b , and c=λ −1/2 c .For each eigenvalue, a eigenvector q can be calculated by the standard procedure of solving the eigenvalue equation Eq=λq, which gives the direction of the corresponding axis.

Fig. 1 .
Figure1shows 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/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 disc-like

Fig. 2 .
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 .
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 .
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 tot in the simulation.

Fig. 5 .
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/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).
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).(C.3) 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 det(E − λI 3 ) = 0, (C.14)
The individual entries are