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.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
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.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
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.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
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.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
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).


Figure 1: A gentle contact between two rubble-pile bodies formerly in a close, synchronous orbit. The shape of the progenitor bodies are preserved, and the neck-like contact area remains narrow and well-defined.               

The results of our suite of simulations suggest that any merger that could produce an object like Arrokoth would have to have been quite slow (roughly equal to the mutual escape speed of the progenitor objects or slower) and at a near-grazing angle. We find that the most plausible scenario is a gradual inspiral of two bodies in a close synchronous orbit leading to final, gentle merger.

Acknowledgements: The simulations described here were performed on the Deepthought 2 cluster at the University of Maryland, College Park. This work was supported by NASA grant NNX15AH90G awarded by the Solar System Workings program, by a University of Maryland Graduate School Research and Scholarship Award, and by NASA's New Horizons project via contracts NASW-02008 and NAS5-97271/TaskOrder30.


[1]: Stern, S.A., Weaver, H.A., Spencer, J.R., Olkin, C.B., Gladstone, G.R., Grundy, W.M., Moore, J.M., Cruikshank, D.P., Elliott, H.A., McKinnon, W.B. and Parker, J.W., 2019. Initial results from the New Horizons exploration of 2014 MU69, a small Kuiper Belt object. Science, 364(6441).

[2]: Marohnic, J.C., Richardson, D.C., McKinnon, W.B., Agrusa, H.F., DeMartini, J.V., Cheng, A.F., Stern, S.A., Olkin, C.B., Weaver, H.A. and Spencer, J.R., 2020. Constraining the final merger of contact binary (486958) Arrokoth with soft-sphere discrete element simulations. Icarus, p.113824.

[3]: Stadel, J.G., 2001. Cosmological N-body simulations and their analysis (Doctoral dissertation, University of Washington).

How to cite: Marohnic, J. C., Richardson, D. C., McKinnon, W. B., Agrusa, H. F., DeMartini, J. V., Cheng, A. F., Stern, S. A., Olkin, C. B., Weaver, H. A., and Spencer, J. R.: Constraining the final merger of contact binary (486958) Arrokoth with soft-sphere discrete element simulations, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-378,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Aldo Dell'Oro, Jacopo Boccenti, Federica Spoto, Paolo Paolicchi, and Zoran Knezevic

The availability of  large data set of asteroids  proper elements made
possible in recent years a deeper and deeper analysis of the structure
of the  asteroid dynamical families in  the Main Belt (Milani  et al.,
2014). At the same time the  Yarkovsky effect, the importance of which
has been  largely recognized,  has imposed  a new  data interpretation
paradigm.   Precisely  the  Yarkovsky  effect  is  the  basis  of  one
estimation method  of the age  of the asteroid families  analyzing the
shape of the family members  distribution in the plane diameter versus
proper semimajor axes, the so called V-shape (Spoto et al., 2015).

This very important  goal has been achieved under  the assumption that
the  post-formation  evolution of  the  semimajor  axes producing  the
V-shape  structure was  linear  in time  and  dominates their  present
dispersion.  In  particular the  impact of the  initial values  of the
fragments ejection  velocities, determining the initial  dispersion of
the semimajor  axes, is supposed  to be marginal,  assumption probably
correct for old enough families. If  this is not the case, the initial
ejection velocity  field introduces a more  or less large bias  in the
age  determination.  We  conducted a  series of  numerical experiments
simulating  the  formation and  dynamical  evolution  of the  asteroid
families in order to evaluate the extent of this bias.

Moreover in  our model we  included a series of  mechanisms concerning
the initial  break-up outcomes.  In particular  two aspects  have been
taken  into   account:  the   effect  of  the   initial  gravitational
re-accumulation  on the  structure of  the V-shape  and the  effect of
possible  correlations between  spin direction  and ejection  velocity
direction (Milani et  al., 2019). We show how the  former effect could
be responsible of  some systematic features of  the V-shape structure,
while  the second  can in  some  cases mitigate  the bias  in the  age
determination introduced by the initial ejection field.

Finally we discuss the role of  the family members spin axes evolution
due to  non destructive  collisions, affecting  the efficiency  of the
Yarkovsky-driven  semimajor axis  mobility.  We show  how this  latter
mechanism modifies the determination of the age.



Milani A., Cellino A., Knezevic Z., Novakovic B., Spoto F., Paolicchi P.
Asteroid families classification: Exploiting very large datasets.
Icarus, 239, 46-73.

Spoto F., Milani A., Knezevic Z.
Asteroid family ages.
Icarus, 257, 275-289.

Milani A., Knezevic Z., Spoto F., Paolicchi P.
Asteroid cratering families: recognition and collisional interpretation.
Astronomy and Astrophysics, 622, A47.

How to cite: Dell'Oro, A., Boccenti, J., Spoto, F., Paolicchi, P., and Knezevic, Z.: Asteroid families age determination: the role of the physical effects, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-90,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Valerio Carruba

Asteroid families are groups of asteroids that are the product of collisions or of the rotational fission of a parent object.  These groups are mainly identified in proper elements or frequencies domains.   Because of robotic telescope surveys, the number of known asteroids has increased from about 10,000 in the early 90's to more than 750,000 nowadays. Traditional approaches for identifying new members of asteroid families, like the hierarchical clustering method (HCM), may   struggle to keep up with the growing rate of new discoveries. Here we used machine learning classification algorithms to identify new family members based on the orbital distribution in proper (a,e,sin(i)) of previously known family constituents. We compared the outcome of nine classification algorithms from stand alone and ensemble approaches.  The Extremely Randomized Trees (ExtraTree) method had the highest precision, enabling to  retrieve up to 97% of family members identified with standard HCM.

How to cite: Carruba, V.: Machine learning classification of new asteroid families members, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-36,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Pierre Deram, Agnès Fienga, and Mickaël Gastineau

Since December 2013, the GAIA spacecraft (ESA) is observing the sky with an unprecedented accuracy. Gaia DR2, released in April 2018, contains the position and epoch of 14099 known solar system objects (SSOs) representing more than 2 million observations collected during the first 22 months of operation. In this presentation, we used the new released INPOP19a planetary ephemerides to perform their orbital adjustment and compare them to radar and optical ground-based observations. In order to reduce the time of computation and in anticipation of the huge amount of datas expected with future DR3, a specific solving strategy for the normal equations of the Gauss-Newton algorithm is presented, making the best use of the model design. We obtain post-fit residuals that are closed to the expected performance of GAIA and overall consistent with the values announced by the DR2 reference orbit determination (see Gaia Collaboration et al 2018). In order to disccus the reliability of the obtained orbits, a combination with radar and optical ground based observations was performed for 23 objects using two different numerical methods: a systematic exploration of the weighting scheme coupled with a residual post-fit analysis, and a Least Square Variance Component Estimation adjustment algorithm (LSVCE). Such methods can be extended to all inverse problems within the framework of least-square formalism.

How to cite: Deram, P., Fienga, A., and Gastineau, M.: Asteroids orbit determination using INPOP19a adjusted on GAIA observations: strategy and validation. , Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-238,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Robert Melikyan, Beth Clark, Carl Hergenrother, Steve Chesley, Michael Nolan, Quanzhi Ye, and Dante Lauretta

NASA’s OSIRIS-REx mission set out in September of 2016 to survey and sample the near-Earth asteroid (101955) Bennu. During observations of Bennu, mission scientists observed centimeter-scale pebbles ejecting off the surface of the asteroid [1]. Many of these particles have been observed to follow hyperbolic trajectories, and the intersection of the orbits of Bennu and Earth suggests the possibility of particle flux at Earth [2]. We simulate the evolution of the motion of particles ejected from Bennu with a focus on potential meteor activity at Earth.


We developed a complex simulated environment that accounts for the most potent perturbing gravitational bodies and solar radiation forces that are applicable to our centimeter-scale size range [3]. We use REBOUND, an orbital integration API developed by Hanno, Rien and Tamayo [4], augmented with REBOUNDx, developed by Tamayo et al. [5], to include solar radiation pressure and Poynting-Robertson (PR) drag. Bennu is modeled as a massless object, though great care is spent ensuring Bennu’s orbital accuracy as the particles are integrated towards well-defined close approaches in the near future [6]. Contrary to most meteoroid stream evolution studies, our simulations release particles from Bennu at a regular cadence (600 grams per week) throughout its orbit to resemble mission observations.


To test the accuracy of the model and integrator, we prepared test cases for expected behavior. Recreating examples of particle resonance traps due to PR drag confirmed the implementation of the non-gravitational forces. Observations of the particle stream circularization (Fig. 1) and associability over time were also indicative of expected behavior. Additionally, planetary bodies were initialized from a JPL ephemeris (DE 431) and continuously compared to these values as a reference. The modeled solar system maintained an accuracy of 10-1 to 101 arcseconds of mean anomaly for their expected positions, giving us confidence in our methods.


Towards the main objective of observing Earth-particle interactions, limitations in processing power and time encourage us to run simulations at particle production rates lower than the mass loss observed at Bennu by a factor of 100. Earth-particle close approaches are recorded annually during simulations. These data are later converted into stream density estimates and Earth impact probabilities through analysis on the B-plane [7]. This formalism enables the calculation of the statistical likelihood of impact with Earth for all of the particles produced in the simulation over the 348 years of interest (1788–2135); the years over which Bennu’s position is best constrained [6]. Such results now inform our conclusions on zenith hourly rate flux measurements at Earth.


This work will be of practical importance for professional and amateur astronomers searching for Bennuid meteors. While the exact particle production mechanisms are still open to debate, we are hopeful that this initial work can be generalized to encompass the entire near-Earth asteroid population. This work, together with additional particle ejection observations and analyses from the OSIRIS-REx mission, will pave the way to a full understanding of this astronomical phenomenon.


[1] D.S. Lauretta and C.W. Hergenrother et al. (2019). Science 366, 1217-1227.

[2] Q. Ye (2019) Notes of the AAS 3, 56.

[3] P. Jenniskens et al. (2011) Icarus 216, 40-61, &

[4] H. Rien and S.F. Liu (2012) Astronomy & Astrophysics 537, A128.

[5] Tamayo, Daniel, et al. (2020) Monthly Notices of the Royal Astronomical Society 491.2

[6] S.R. Chesley et al. (2014) Icarus 235, 5-22

[7] D. Farnocchia et al. (2019) Springer, Celestial Mechanics and Dynamical Astronomy.

How to cite: Melikyan, R., Clark, B., Hergenrother, C., Chesley, S., Nolan, M., Ye, Q., and Lauretta, D.: A Natural Sample Delivery Mechanism: Estimating the Flux of Bennuid Meteors, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-428,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Stavro L. Ivanovski, Alice Lucchetti, Maurizio Pajola, Ivano Bertini, Giovanni Zanotti, Davide Perna, Elisabetta Dotto, Vincenzo Della Corte, Marilena Amoroso, Simone Pirrotta, Andrea Capannolo, Michele Lavagna, Alessandro Rossi, Eugene G. Fahnestock, Masatoshi Hirabayashi, Sabina D. Raducan, Andrea Meneghin, Giovanni Poggiali, John R. Brucato, and Gabriele Cremonese and the and the LICIACube team


On September 30, 2022, NASA’s Double Asteroid Redirection Test (DART, [1]) mission will be the first space mission demonstrating the kinetic impactor method for planetary defense. DART will impact Didymos B (a 164±18 m secondary of asteroid (65803) Didymos [2]) with its mass of 650 kg and at a speed of 6.6 km/s. Such impact is expected to change the secondary orbital period by about 10 minutes. DART will carry as a piggyback the Light Italian Cubesat for Imaging of Asteroids (LICIACube, [3]) which will be released from DART ten days before the impact. LICIACube will provide evidence of the impact and will take multiple images of the target up to a distance of ~55 km from the target. The LICIACube narrow and wide angle cameras - LEIA (LICIACube Explorer Imaging for Asteroid) and LUKE (LICIACube Unit Key Explorer), respectively – will capture the post-impact processes coming from in situ events, such as the newly formed crater, the expanding ejecta and the dynamics of its plume.  In particular, the measurement of the motion of the slow (<5 m/s) ejecta in the plume will be feasible with LICIACube, thus allowing description of its structure and the evolution of the dust size and velocity distribution.  In this work, we study the motion of the plume evolution in terms of dust numerical simulations. To accomplish this, we modified the non-spherical dust model [4] and applied it to an asteroid impact scenario using analytical expressions to estimate the initial velocity and mass of ejecta.  As input we used the ejecta impact properties (ejecta mass, velocity, launch position distribution, orientation) constrained with iSALE numerical simulations [5]. We discuss the influence of the non-sphericity of the particles on the dynamical properties of the plume, such as the velocity and dust spatial distribution, and address the optical thickness not only in terms of particle size distribution but also as a function of particle shape and orientation.

Dust plume model

We use a 3D+t non-spherical dust model that solves the Euler dynamical and kinetic equations. Considering free-collisional dust regime we study the effects on the particle dynamics provided with different shapes, initial particle orientation and velocities as well as torque. Torque is computed from the law of variation of the angular momentum by using the Euler dynamic equations. The particles are assumed to be homogeneous, isothermal convex bodies, having the same physical properties of the target. The dust motion is governed by solar radiation pressure force, initial dynamical parameters (speed, orientation and torque) and gravity of the asteroid binary system. To compute the gravity we take into account that Didymos B is orbiting at distance of ~1.18 km the Didymos primary (780±80 m in size) [6].

As input we use the physical parameters obtained via numerically simulated impact into low-gravity and strength- dominated asteroid surface done with the iSALE numerical code [5]. A detailed description of the parameters set used in the impact modelling are reported in Table 1. As an example, in Fig.1, we show the iSALE simulation of the spatial distribution of tracer particles for 50 ms after the impact.

Table 1: Paremeters used for iSAle simulation. A detailed description can be found in [5].

Fig. 1: A snapshot of the formation of the crater resulting  from the iSALE modelling performed with parameters reported in Table 1. The colorbar represent the peak pressure of the ejecta tracer particles.

Non-spherical dust motion in the plume

We have studied the motion of non-spherical particles after the impact assuming various set of initial parameters coming from the iSALE simulations, such as particle mass, shape, initial orientation, launch position and computed velocity. The Didymos system physical parameters and particle temperature were taken from [7]. Here, in Fig. 2, we show the velocity and the distance from the surface where it has been reached for a prolate spheroid (aspect ratio = 2 [4]) with assumed density 2600 kg/m3. The particle mass is 7.45 x10-2 kg, the initial velocity 0.94 m/s, the launch position is 11 m from the center of the crater and the initial tilt of the initial velocity vector away from local surface normal is 46 degrees. We plotted the results for two cases: 1) only the Didymos B gravity field is considered and 2) the gravity of the whole binary is used. Our results show that the particle dynamics is sensitive to parameters such as orientation and ejection velocity.  

Fig. 2: Velocity and distance at which it has been achieved for a prolate particle (aspect ratio a/b = 2). The complete set of parameters is given in the text. The solid line indicates the case when only gravity of Didymos B  was considered and the dashed line when the gravity of the whole binary was used (binary was approximated as a mass point).

This study of dust dynamics in expanding ejecta together with LICIACube plume evolution images can further constrain the momentum transfer efficiency estimation as well as the ejecta velocity and dust size distribution in the “far-field” expanding plume.



This research was supported by the Italian Space Agency (ASI) within the LICIACube project (ASI-INAF agreement AC n. 2019-31-HH.0).



 [1] Cheng A. F. et al. (2016) Planet. Space Sci., 121, 27-35. [2] de León, J. et al. (2010) Astron. Astrophys, 517, A23. [3] Dotto, E. et al., (submitted), Planet. Space Sci. [4] Ivanovski S. et al. (2017), Icarus 282, p. 333 - 350. [5] Raducan, S. D. et al. (2019), Icarus, 329, p. 282 - 295. [6] Michel, P. et al. (2016) Adv. Space Sci. 57, p. 2529 - 2547. [7] Yu Y. et al. (2017), Icarus 282, p.313 - 325

How to cite: Ivanovski, S. L., Lucchetti, A., Pajola, M., Bertini, I., Zanotti, G., Perna, D., Dotto, E., Della Corte, V., Amoroso, M., Pirrotta, S., Capannolo, A., Lavagna, M., Rossi, A., Fahnestock, E. G., Hirabayashi, M., Raducan, S. D., Meneghin, A., Poggiali, G., Brucato, J. R., and Cremonese, G. and the and the LICIACube team: Modelling dust distribution in the ejecta plume from nonspherical dust dynamics perspectives in support of the LICIACube and DART missions, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-1096,, 2020.

Chairperson: Stavro Ivanovski, Vladimir Zakharov
Dorothea Bischoff, Bastian Gundlach, and Jürgen Blum

In the Solar System, several different types of small bodies can be found, whose histories can deviate highly from each other, including time and process of formation, followed by collisional or other evolutionary processes. For cometary nuclei, the history and formation remain unsolved and different hypotheses exist, which may explain their past (Weissman et al. 2020). In this study here, the focus lies on the thermophysical properties of small bodies. Different formation scenarios lead to different properties of the resulting bodies. In the scenario of the gravoturbulent collapse of a pebble cloud, macroscopic void spaces remain in the body (e.g. Blum et al. 2014, 2017). In contrast, the hierarchical accretion leads to more compact bodies without large voids (e.g. Davidsson et al. 2016). This effects the thermophysical behaviour, as in the first case, besides the network conductivity of the solid material, also radiation can efficiently transport heat through the voids. The radiative heat conductivity can be described as follows:

λrad = 16/3 σ TΛ(Φ)

where σ, T, Λ and Φ are the Stefan-Boltzmann constant, the temperature, the radiative mean free path and the volume filling factor, respectively. The network part of the heat conductivity can be written as

λnet = H λagg/par(T)


with the Hertz factor H and the aggregate or particle thermal conductivity λagg/par (see Gundlach and Blum 2012, and Gundlach et al 2020 for details). The radiative thermal conductivity has a much higher dependency on temperature than the network part.

Here, we address the question, how radiation influences the thermal behaviour of a planetesimal in a measurable way. Therefore, the following steps are foreseen: a thermophysical model is set up, which calculates the heat transport in a small body that is illuminated by the Sun. For simplicity, we assume a circular orbit around the Sun at several heliocentric distances without obliquity and varying rotational periods. The model solves the heat transport equation in a one-dimensional approximation. In this study, we focus on the surface temperature as the output. This is motivated by the measurement capabilities of most spacecraft visiting small bodies. To measure the surface temperature is relatively easy, often done and enough accurate compared to the remote measurements of inner temperatures.

Figure 1 : Results of the thermophysical modelling for 2 AU (left) and 40 AU (right): surface temperature over time for the simulation with (dashed) and without (continuous) radiative part.

With the described model, the influence of the radiative part of the thermal conductivity is studied by varying it in different setups, namely heliocentric distance, and rotation period. Thereby, the Fourier number, which represents the numerical accuracy of the model run, is optimized to be near its favourable value of 0.5. In a first attempt, only the surface temperature of the equator was regarded. An exemplary plot of two model runs is shown in Figure 1 for a rotation period of 10 hours and at heliocentric distances of 2 AU and 40 AU, respectively. For each heliocentric distance, the maximal temperature is similar for the two cases (with and without radiative conductivity) studied, but during the night, cooling behaves differently. In the case of 2 AU, a median thermal conductivity of 0.02 W/(m K) was calculated; for 40 AU, the median thermal conductivity reduces to 0.003 W/(m K). In this colder regime, the differences at night are reduced. In a next step, this analysis will be expanded to the surface temperature of the meridian from pole to equator, assuming a spherical body.

The first results indicate that for some parameter sets, the two scenarios show sufficient differences to be distinguishable by measurement. Comparing the results of our study to measured data from small bodies could allow further insights into the formation or evolutionary paths of small bodies. More detailed results and conclusions will be presented at the EPSC 2020.



This work was funded through the DFG project BL 298/27-1.



Blum J., Gundlach B., Mühle S., Trigo-Rodriguez J. M., 2014, Icarus, 235, 156

Blum J. et al., 2017, MNRAS, 459, S755

Davidsson B. J. R. et al., 2016, A&A, 592, A63

Gundlach B., Blum J., 2012, Icarus, 219, 618

Weissman P., Morbidelli A., Davidsson B., Blum J., 2019, Space Sci. Rev., 216, 6

How to cite: Bischoff, D., Gundlach, B., and Blum, J.: Thermophysical modelling of small bodies: influence of the radiative thermal conductivity, Europlanet Science Congress 2020, online, 21 Sep–9 Oct 2020, EPSC2020-951,, 2020.