Europlanet Science Congress 2020
Virtual meeting
21 September – 9 October 2020
Europlanet Science Congress 2020
Virtual meeting
21 September – 9 October 2020

Oral presentations and abstracts


The goal of this session is to cover numerical simulations and relevant laboratory investigations related to the Small Bodies (comets, KBOs, rings, asteroids, meteorites, dust), their formation and evolution, and the instruments of their exploration. This session is specially focused on the interdisciplinary approach in the development of models (formal descriptions of physical phenomena), experiments (on ground and in micro-gravity), and mathematical simulations (computational methods and algorithms of solution) of various astrophysical phenomena: (i) dusty gas cometary atmospheres; (ii) volcanic activity on icy satellites (e.g. Enceladus and Io); (iii) planetary body formation (e.g. via pebbles growth), and planetesimal dynamics.

This session will include an introduction and discussion of new and/or existing laboratory studies in simulated space-like environments and models, experimental techniques, computational methods that can address the results of analytical,experimental and numerical analysis (with respect to computational methods and algorithms of solution) on the above described studies.
Abstracts on thermophysical evolution models of small bodies interiors as well as on the modeling of atmosphere and exosphere are welcome.

Conveners: Stavro Lambrov Ivanovski, Vladimir Zakharov | Co-conveners: Vincenzo Della Corte, Michelangelo Formisano, Marco Fulle, Raphael Marschall, Luis Diego Pinto, Alessandra Rotundi, Diego Turrini

Session assets

Session summary

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Munan Gong, Alexei Ivlev, Bo Zhao, and Paola Caselli

Grain growth in protoplanetary disks is the first step towards planet formation. One of the most important pieces in the grain growth model is calculating the collisional velocity between two grains induced by the turbulent motion of the gas. The collisional velocities in previous works are obtained based on the assumption that the turbulence is hydrodynamic with the Kolmogorov power spectrum. However, realistic protoplanetary disks are magnetized, and turbulent motions can be induced by the magneto-rotational instabilities (MRI).

To understand the impact of MRI turbulence on grain growth, we carry out 3D MHD simulations of the MRI as well as driven turbulence, for a range of physical and numerical parameters. Figure 1 illustrates the magnetic field structure in one of the MRI simulations. We investigate two fundamental parameters of the turbulence that can have large impacts on grain collision velocities, the kinetic energy spectrum and the turbulence autocorrelation time. 

We find that the convergence of the turbulence α-parameter does not necessarily imply the convergence of the energy spectrum. The MRI turbulence is largely solenoidal, for which we observe a persistent kinetic energy spectrum with a -4/3 slope, shallower than the -5/3 slope in the Kolmogorov turbulence (Figure 2). The same is obtained for solenoidal driven turbulence with and without magnetic field, over more than 1 dex near the dissipation scale. This power-law slope appears to be converged in terms of numerical resolution, and to be due to the bottleneck effect. The kinetic energy in the MRI turbulence peaks at the fastest growing mode of the MRI. In contrast, the magnetic energy peaks at the dissipation scale. The magnetic energy spectrum in the MRI turbulence does not show a clear power-law range, and is almost constant over approximately 1 dex near the dissipation scale. The turbulence autocorrelation time is nearly constant at large scales, limited by the shearing timescale, and shows a power-law drop at small scales, with a slope close to -1, steeper than that of the eddy crossing time.

The deviation from the standard picture of the Kolmogorov turbulence can have a significant impact on the grain collisional velocities. Figure 3 shows the collisional velocities in Kolmogorov turbulence and turbulence observed in the simulations. For small grains that are well coupled to the gas, the collisional velocities can vary by orders of magnitude across different turbulence models. To understand the dust growth in protoplanetary disks, better computational techniques are needed to achieve the large dynamical range necessary for resolving the turbulence spectrum in the future.


Gong M., Ivlev A. V., Zhao B., Caselli P., 2020, ApJ, 891, 172

How to cite: Gong, M., Ivlev, A., Zhao, B., and Caselli, P.: Impact of magneto-rotational instability on grain growth in protoplanetary disks, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-316,, 2020.

James E. Robinson, Wesley C. Fraser, Alan Fitzsimmons, and Pedro Lacerda


In the trans-Neptunian region a relatively large fraction of minor bodies are found to exist as mutually bound binary pairs. The binary fraction varies for different dynamical classes and is highest for the Cold Classicals, which are dynamically unexcited planetesimals generally found at 39 < a < 48 AU, with estimates of up to 30% of the Cold Classicals being binaries (Noll et al. 2008).

A distinct feature of the observed trans-Neptunian binaries is that they often have components that are of a similar size. These components also have similar optical colours which likely indicates similar surface compositions. Furthermore, the binary orbits span a range of separations but it is not uncommon for them to be found on wide (abin > 0.05 RHill) dynamically fragile orbits. The binary orbit inclinations range from prograde to retrograde, with ~80% of the known trans-Neptunian binaries having prograde orbits. The discovery of the bi-lobate nature of the Cold Classical contact binary 2014 MU69 Arrokoth, with its similar size and colour lobes, indicates that this object is likely the end state of a previously separated trans-Neptunian binary (Stern et al. 2019).

Previous work (Nesvorný et al. 2010) has shown that the formation of planetesimals via gravitational collapse is a promising mechanism that can reproduce the unique properties of the trans-Neptunian binaries. We assume that a slowly rotating cloud of pebbles grows in the protoplanetary disk, for example via the streaming instability. Eventually this cloud collapses under its own self gravity and merging collisions lead to the rapid growth of ~100 km planetesimals. Conservation of angular momentum in the rotating cloud means that these planetesimals are frequently formed as a binary pair. For a given binary orbit the angular momentum is maximised by forming a system with similar mass components. Additional work (Nesvorný et al. 2019) has shown that the inclination distribution of binary planetesimals formed via streaming instability and gravitational collapse is a good match for the observed inclination distribution.


We performed N-body simulations of the gravitational collapse of a pebble cloud using the Symplectic Epicyclic Integrator within the REBOUND package (Rein & Liu 2012). The pebble cloud was approximated by a rotating uniform spherical distribution of 105 computational particles. The radii of the computational particles was inflated in order to replicate the collisional behaviour of a high number density pebble cloud. Collisions were detected when particles were found to be overlapping and we considered the timestep required to ensure that collision detection was robust. We performed a deep search for systems of gravitationally bound particles and tested their stability. For all systems produced we investigated the population characteristics of their mass and orbital parameters.


Our results demonstrate that gravitational collapse is an efficient producer of bound planetesimal systems, and frequently forms multiple bound systems per cloud. On average there were ~1.5 bound systems produced per cloud in the mass range we have studied. As well as the large equal-sized binaries which were the focus of previous work, we found that gravitational collapse can produce a range of systems such as massive bodies with small satellites and low mass binaries with a high mass ratio. The binaries with equal-size components have bound orbits with low inclinations and low to moderate eccentricity, similar to the observed trans-Neptunian binaries. The additional binaries with dissimilar-sized components or low system mass span a large range of inclinations and eccentricities. Our results disfavour the collapse of the high mass clouds in our dataset (Mc = 1.8e21 kg), as they form large equal-sized binaries that are not observed (figure 1). This finding is in line with reported upper mass limits of clouds formed by the streaming instability. Gravitational collapse of low mass clouds (Mc ≤ 4.2e18 kg) can create binary systems with similar total mass and mass ratio as the contact binary Arrokoth (figure 1). Furthermore, we find that collisions between planetesimals growing in such a collapsing cloud should be gentle enough to preserve a bi-lobed structure, which further supports the gravitational collapse origin of Arrokoth.


Figure 1. Primary V band magnitude against magnitude difference between primary and secondary components for binaries produced by gravitational collapse simulations. Marker colour and shape indicates the initial cloud mass. The vertical line indicates the mass ratio cut off of m2/m1 ≥ 10−3 in the orbit search algorithm. The observed trans-Neptunian binaries from Grundy (2019) are shown as black 'x's. We represent 'special cases' and dwarf planets as red '+'s. The red circle with error-bars represents how Arrokoth would appear if its components could be separately resolved. An approximate empirical detection limit (Noll et al. 2008) and a lower magnitude limit of 25 are shown as red dotted lines.


Grundy, W. (2019). Mutual Orbits of Binary Transneptunian Objects. Retrieved December 2nd, 2019, from online webpage website:

Nesvorný, D., Youdin, A. N., & Richardson, D. C. (2010). Formation of Kuiper Belt Binaries by Gravitational Collapse. The Astronomical Journal, 140, 785–793.

Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. (2019). Trans-Neptunian binaries as evidence for planetesimal formation by the streaming instability. Nature Astronomy, 3, 808–812.

Noll, K. S., Grundy, W. M., Chiang, E. I., Margot, J.-L., & Kern, S. D. (2008). Binaries in the Kuiper Belt. In The Solar System Beyond Neptune (pp. 345–363). Tucson: The University of Arizona Press.

Rein, H., & Liu, S. F. (2012). REBOUND: An open-source multi-purpose N-body code for collisional dynamics. Astronomy and Astrophysics, 537, A128.

Stern, S. A., Weaver, H. A., Spencer, J. R., Olkin, C. B., Gladstone, G. R., Grundy, W. M., … Zurbuchen, T. H. (2019). Initial results from the New Horizons exploration of 2014 MU69, a small Kuiper Belt object. Science, 364(6441).







How to cite: Robinson, J. E., Fraser, W. C., Fitzsimmons, A., and Lacerda, P.: Investigating Gravitational Collapse of Pebble Clouds to form Trans-Neptunian Binaries, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-891,, 2020.

Sergei Ipatov

   Migration of planetesimals to the Earth from the zone beyond the orbit of Jupiter was considered by us e.g. in [1-3]. References to papers of several other authors on migration of bodies to the Earth were presented in [1]. In [1] we considered migration of planetesimals from the zone from 4.5 to 12 AU. In [2-3] migration of planetesimals with initial semi-major axes a0 of their orbits between 5 and 40 AU was considered. Below I also consider migration of planetesimals with a0 between 3 and 5 AU to the Earth.

   Migration of planetesimals under the gravitational influence of 7 planets (from Venus to Neptune) or 5 planets (from Venus to Saturn) was calculated with the use of the symplectic code from [4]. The. In each variant of the calculations, the initial values of semimajor axes of orbits of planetesimals varied from amin to amax=amin+da, the initial eccentricities were equal to eo, and the initial inclinations equaled to eo/2 rad. Orbital elements of the migrated planetesimals were recorded in computer memory with a step of 500 years. Based on these arrays of orbital elements, I calculated the probabilities of collisions of planetesimals with the Earth, and for some runs I also calculated the probabilities of collisions of the planetesimals with other terrestrial planets, the Moon and their embryos. The calculations were made similar to those in [1-3, 5-7].

    In the series of calculations considered in [2-3], da=2.5 AU, and amin took values from 2.5 to 40 AU in increments of 2.5 AU. The initial eccentricities equaled to 0.3 or 0.05. In each calculation variant, 250 planetesimals were considered, but for the same values of amin, da, and eo, several (up to 8) calculation variants were performed. So the total number of considered planetesimals for a set with fixed values of amin, da and e0 could reach 2000. Some calculations were made for da=0. In the recent series of calculations, da=0.1 AU, and amin took values from 3.0 to 4.9 AU in increments of 0.1 AU. For this series of calculations, the initial eccentricities equaled to 0.02 or 0.15. In Figs. 1-6 for several series of calculations, we present the values of 106pE, where pE is the probability of a collision of a planetesimal with the Earth. In Figs. 1-2 the values of 106pE are presented for time intervals equaled to 1, 10 and 100 Myrs. For other figures usually greater time intervals (up to 2 Gyrs) were considered (until the values of pE finished to grow with time). For variants presented in Figs. 1-4 and 6, the gravitational influence of 7 planets (from Venus to Neptune) was taken into account. For runs for Fig. 5, Uranus and Neptune were excluded.

   At amin≤10 AU, the value of pE calculated for a run with 250 bodies could vary hundreds of times for different calculation variants with the same values of amin, da and eo. Such difference was earlier found for calculations of migration of Jupiter-crossing objects [5-6]. One among several hundreds or among thousands of such objects moved in Earth-crossing orbits during millions or even tens of millions of years, and the probability of a collision of such object with the Earth was greater than that for hundreds or even thousands of other objects. The values of pE in Figs. 1-6 vary from a value less than 10-7 to values of the order of 10-3, but they are mainly between 10-6 and 10-5. On average, pE is smaller for greater amin. There were no runs with pE>10-5 at amin≥12.5 AU or for the runs without Uranus and Neptune. pE≤1.5×10-6 at e0=0.3 and amin≥27.5 AU. In Figs. 1-2, the fraction of runs with pE>10-5 at e0=0.15 is greater than at e0=0.02. In some runs there was a considerable growth of pE after 10 Myr. At da=0.1 AU and e0=0 or e0=0.15, there were more runs with pE>2×10-5 for 3.2≤amin≤4.1 AU than for other values of amin. 

   Calculations showed that the amount of material delivered from beyond Jupiter’s orbit to the Earth could exceed the mass of the Earth’s oceans if the total mass of planetesimals beyond Jupiter’s orbit was about 200 Earth masses. Some (perhaps 1/3) of this material consisted of water and volatile substances. The mass of the substance delivered to the planet to the mass of the planet for Mars was approximately two times greater than for Earth, and such relations for Mercury and Venus were slightly larger than for Earth. The total mass of planetesimals migrating from beyond the orbit of Jupiter and colliding with the Moon was 16 or 17 times less than the total mass such bodies collided with the Earth.

   The work was carried out as a part of the state assignments of the Vernadsky Institute of RAS.


   [1] Marov M.Ya., Ipatov S.I., Delivery of water and volatiles to the terrestrial planets and the Moon // Solar System Research, 2018, 52, 392-400.

   [2] Ipatov S.I., Migration of bodies to the Earth and the Moon from different distances from the Sun // The Ninth Moscow Solar System Symposium 9M-S3, 2018,, 9MS3-SB-11, p. 104-106,

   [3] Ipatov S.I. Migration of planetesimals to the Earth and the Moon from different distances from the Sun // 50th LPSC, 2019, #2594. . e-poster:

   [4] Levison H.F., Duncan M.J. The long-term dynamical behavior of short-period comets // Icarus, 1994, 108, 18-36.

   [5] Ipatov S.I., Mather J.C. Migration of Jupiter-family comets and resonant asteroids to near-Earth space // Annals of the New York Academy of Sciences, 2004, 1017, 46-65.  

   [6] Ipatov S.I., Mather J.C. Comet and asteroid hazard to the terrestrial planets // Advances in Space Research, 2004, 33, 1524-1533.

   [7] Ipatov S.I., Probabilities of collisions of planetesimals from different regions of the feeding zone of the terrestrial planets with the forming planets and the Moon // Solar System Research, 2019, 53, 332-361.




How to cite: Ipatov, S.: Migration of planetesimals from beyond Mars’ orbit to the Earth, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-71,, 2020.

Joseph DeMartini and Derek Richardson

It is believed that many near-Earth asteroids are rubble piles: aggregates of smaller regolith particles bound together by self gravity and cohesive forces, formed from the reaccumulation of much larger, monolithic progenitors after catastrophic collisions [1]. The recent OSIRIS-REx and Hayabusa2 sample-return missions to Bennu and Ryugu, respectively [2, 3], have shown complex and diverse regolith surfaces on rubble piles, with grains spanning from millimeters to tens of meters in radius, evidence of granular flow (slumping), and artifacts such as large surface boulders without surrounding craters (potentially indicating vertical regolith migration).

A successful approach to simulating rubble piles is the discrete element method (DEM), which uses individual (typically spherical) particles to represent the constituent blocks of the rubble pile, with particle properties (including friction parameters) chosen to mimic bulk properties of real materials, such as angle of repose. Particles that make up real granular materials, however, are most often not spherical. The shapes of real grains make it more difficult for them to slide past one another, can create quasi-stable equilibria when grains are vertically stacked, and cause the material to “bulk” (creating more void space and less efficient packing). These phenomena of real granular materials are not well exhibited by spherical particles, thus the most accurate way to capture a realistic picture of granular dynamics is to use irregularly shaped particles instead of spheres in simulations. The aim of this work is to present a new algorithm for handling simulations using large numbers of irregularly shaped particles with pkdgrav, a DEM code that has been used effectively in the past to model granular dynamics on rubble piles [4, 5, 6, 7], and show that simulating with irregularly shaped grains can lead to results more akin to real granular systems.

We construct irregular particles by “gluing” spheres together to make bonded (rigid) aggregates. Originally, aggregate routines for pkdgrav were written to handle a few large aggregates made up of many individual particles. In this scheme, aggregate properties are tracked in serial, where a brute force search finds particles belonging to an aggregate then operates on them, searching through every particle on every processor for every aggregate in the simulation. An inefficiency arises when the number of aggregates is similar to the total number of particles (N), as in a simulation of a granular bed made up of small aggregates, each only containing a few particles. Here, the brute force search becomes an ∼O(N2) operation. This is prohibitively expensive for processes requiring large N, like granular dynamics simulations.

We solved the issue of efficiently locating particles that belong to an aggregate by forcing a reordering of the particles before the search ever occurs, to put them in order of their aggregate index numbers on each processor. We then exploit this particle ordering by replacing the original brute-force search with a binary-search algorithm [8] that scales as O(NlogN) and added a “cache-line” method. The binary search is used whenever the particle information is not immediately known, i.e., not in the cache (like the first time a search needs to occur in a simulation). The cache line stores the index of the previously found particle in order to easily step forward and find the next particle that needs to be acted on when the code returns to the function—this is a linear operation in N (O(N) scaling), so it scales well. Preliminary tests show a 25–50% decrease in total (wallclock) runtime for moderate resolution simulations using this algorithm overthe original scheme, and expect better performance at higher resolutions.

To show the importance of being able to simulate irregularly shaped granular material, we apply our method to an investigation of a scientific application for granular dynamics in low gravity: the Brazil-nut effect. The Brazil-nut effect (BNE) [5, 6, 9] is a mechanism for the vertical migration of boulders in a granular medium that has been suggested as an explanation for exposed boulders on asteroid surfaces. In granular dynamics, the BNE is a method in which frictional interactions between particles in a granular medium cause larger blocks to rise to the surface when the medium is subjected to repeated seismic shaking (like Brazil nuts rising to the top of a shaken jar of mixed nuts).

In previous studies with spherical particles, it was shown that increased interparticle friction accounted for a more efficient rise of the large subsurface intruder, but in our preliminary simulations, when simulating with irregular grains, we show that interparticle friction makes the void-filling mechanism less efficient and thus the Brazil nut rises slower. We will explore a broad parameter space in particle size ratio, grain shape, interparticle friction, and interparticle cohesion, to establish how more realistic, irregularly shaped particles exhibit the BNE in low-gravity environments.

This work was supported in part by the GEODES grant, number 80NSSC19M0216, awarded by the NASA SSERVI Program. Simulations were carried out at the University of Maryland on the yorp cluster administered by the Department of Astronomy.
[1] P. Michel, W. Benz, P. Tanga, D. C. Richardson, Science 294 (2001) 1696–1700.
[2] D. S. Lauretta, D. N. DellaGiustina, C. A. Bennett, et al., Nature 568 (2019) 55–60.
[3] S. Watanabe, M. Hirabayashi, N. Hirata, et al., Science 364 (2019) 268–272.2
[4] Y. Yu, D. C. Richardson, P. Michel, S. R. Schwartz, R.-L. Ballouz (2014).arXiv:1408.0168.
[5] S. Matsumura, D. C. Richardson, P. Michel, S. R. Schwartz, R. L. Ballouz, MNRAS 443 (2014) 3368–3380.arXiv:1407.2748.
[6] C. Maurel, R.-L. Ballouz, D. C. Richardson, P. Michel, S. R. Schwartz, MNRAS 464 (2017) 2866–2881.
[7] Y. Zhang, D. C. Richardson, et al. (2017).arXiv:1703.01595.
[8] J. L. Bentley, Communications of the ACM 18 (1975) 509–517.
[9] A. Rosato, K. J. Strandburg, F. Prinz, R. H. Swendsen, Physical Review Letters 58 (1987) 1038–1040.

How to cite: DeMartini, J. and Richardson, D.: An Efficient Algorithm for Operating on Large Numbers of Aggregate Particles with Applications to Simulating the Dynamics of Irregularly Shaped Grains, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-383,, 2020.

Julian C. Marohnic, Derek C. Richardson, William B. McKinnon, Harrison F. Agrusa, Joseph V. DeMartini, Andrew F. Cheng, S. Alan Stern, Cathy B. Olkin, Harold A. Weaver, and John R. Spencer

In 2019, the New Horizons spacecraft flew by the Kuiper belt object (486958) Arrokoth, a member of the cold classical population in the Kuiper belt. Images returned from the New Horizons flyby reveal that Arrokoth is a contact binary composed of two distinct lobes joined by a narrow contact region (the "neck"). The lobes are roughly similar to oblate spheroids in shape and the long axis of the whole body is about 35 km long [1]. Cold classical Kuiper belt objects like Arrokoth are thought be among the most primitive bodies in our solar system, forming in their current positions from the solar nebula over 4 Gyr ago and remaining nearly undisturbed since. The mechanism by which Arrokoth formed is thus of great interest and may lend unique insight into the formation of terrestrial bodies in the outer solar system. Under the assumption that Arrokoth is the union of two distinct progenitor bodies, we use a numerical code to model the final stage of the merger that created this unique object. We investigate a range of possible merger parameters and use the results to determine the most plausible scenarios. We describe the methods we used to simulate the merger and detail our results. The work presented here is described in greater detail in Marohnic et al. 2020, in press for the Icarus special issue on Kuiper belt science [2].

We use the parallel N-body code 'pkdgrav' to model the Arrokoth merger [3]. pkdgrav uses the soft-sphere discrete element method (SSDEM) to compute collisions between spherical particles. SSDEM resolves contacts by allowing particles in contact to interpenetrate slightly as a proxy for deformation. Restoring forces are modeled as damped springs and static, rolling, and twisting friction are accounted for as well. Both lobes are modeled separately, as "rubble piles" composed of spherical particles. In total, each simulated merger event uses between 100,000 and 200,000 particles. Our simulations test the effects of varying impact angle and speed, material cohesion strength, body shape, bulk density, and spin (see Figure 1).